Repository: stfc/PSyclone
Pull request #2962 (2025-04-25 10:35:55)

(Towards #2935) Fix lbound bug when passing array arguments to routines.
This PR started life as per its title but I encountered a section of code (in `sympy_writer`) that needed to be changed to use reference_accesses. In doing that I found there were issues when the same Symbol was accessed in different ways, e.g. `a + a(i)` would result in a crash. Fixing this proved complex for structures but the resulting code is actually not too bad (I claim).


Patch:
diff --git a/changelog b/changelog
index f7892c127..e48884cfd 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,6 @@
+   40) PR #2962 for #2935. Errors in use of LBOUND and looping in generated
+   adjoint kernel.
+
    39) PR #3021. Remove hoisted reference from chunking/tiling
    transformations.
 
diff --git a/doc/developer_guide/sympy.rst b/doc/developer_guide/sympy.rst
index 6d44c048d..b0bc150c4 100644
--- a/doc/developer_guide/sympy.rst
+++ b/doc/developer_guide/sympy.rst
@@ -140,9 +140,12 @@ expression ``a(i:j:k)`` will become
 ``a(i, j, k)``, and to then maintain the same number of arguments for
 each use of an array/function, ``a(i)`` will become ``a(i,i,1)``, and
 ``b(i,j)`` becomes ``b(i,i,1,j,j,1)`` etc. Array expressions like ``a(:)``
-will be using specific names for the lower and upper bound, defaulting
+will use specific names for the lower and upper bound, defaulting
 to ``sympy_lower`` and ``sympy_upper``. So the previous expression
-becomes ``a(sympy_lower, sympy_upper, 1)``. Note that in case of a name
+becomes ``a(sympy_lower, sympy_upper, 1)``. Array accesses that do not have
+parentheses (i.e. just ``a``) are treated similarly but use ``sympy_no_bounds``
+for the lower and upper bounds so ``a`` becomes
+``a(sympy_no_bounds, sympy_no_bounds, 1)``. Note that in case of a name
 clash SymPyWriter will change the names of the boundaries to be unique.
 
 Fortran-specific Syntax
diff --git a/src/psyclone/psyad/adjoint_visitor.py b/src/psyclone/psyad/adjoint_visitor.py
index b44560bd3..ae8f99028 100644
--- a/src/psyclone/psyad/adjoint_visitor.py
+++ b/src/psyclone/psyad/adjoint_visitor.py
@@ -232,11 +232,11 @@ class AdjointVisitor(PSyIRVisitor):
         :param node: an Assignment PSyIR node.
         :type node: :py:class:`psyclone.psyir.nodes.Assignment`
 
-        :returns: a list of PSyIR nodes containing the adjoint \
+        :returns: a list of PSyIR nodes containing the adjoint
             of this node.
         :rtype: list of :py:class:`psyclone.psyir.nodes.Node`
 
-        :raises VisitorError: if the schedule_node method has not been \
+        :raises VisitorError: if the schedule_node method has not been
             called previously.
 
         '''
@@ -247,8 +247,9 @@ class AdjointVisitor(PSyIRVisitor):
                 "as the latter sets up the active variables.")
         assign_trans = AssignmentTrans(self._active_variables)
         new_node = node.copy()
-        # Temporary parent schedule required by the transformation.
-        dummy_schedule = Schedule()
+        # Temporary parent schedule required by the transformation. Needs
+        # access to the ancestor scoping region in order to resolve symbols.
+        dummy_schedule = Schedule(parent=node.scope)
         dummy_schedule.children.append(new_node)
         assign_trans.apply(new_node)
         return dummy_schedule.pop_all_children()
diff --git a/src/psyclone/psyir/backend/fortran.py b/src/psyclone/psyir/backend/fortran.py
index 721d4db44..97f958a28 100644
--- a/src/psyclone/psyir/backend/fortran.py
+++ b/src/psyclone/psyir/backend/fortran.py
@@ -1729,21 +1729,24 @@ class FortranWriter(LanguageWriter):
                 result_list.append(self._visit(child))
         return ", ".join(result_list)
 
-    def call_node(self, node):
+    def call_node(self, node) -> str:
         '''Translate the PSyIR call node to Fortran.
 
         :param node: a Call PSyIR node.
         :type node: :py:class:`psyclone.psyir.nodes.Call`
 
         :returns: the equivalent Fortran code.
-        :rtype: str
 
         '''
         args = self._gen_arguments(node)
-        if isinstance(node, IntrinsicCall) and node.routine.name in [
-                "ALLOCATE", "DEALLOCATE"]:
-            # An allocate/deallocate doesn't have 'call'.
-            return f"{self._nindent}{node.routine.name}({args})\n"
+        if isinstance(node, IntrinsicCall) and node.routine.name not in [
+                "DATE_AND_TIME", "SYSTEM_CLOCK", "MVBITS", "RANDOM_NUMBER",
+                "RANDOM_SEED"]:
+            # Most intrinsics are functions and so don't have 'call'.
+            if not node.parent or isinstance(node.parent, Schedule):
+                return f"{self._nindent}{node.routine.name}({args})\n"
+            return f"{node.routine.name}({args})"
+
         if not node.parent or isinstance(node.parent, Schedule):
             return f"{self._nindent}call {self._visit(node.routine)}({args})\n"
 
diff --git a/src/psyclone/psyir/backend/sympy_writer.py b/src/psyclone/psyir/backend/sympy_writer.py
index 5b261638d..b0fef86fc 100644
--- a/src/psyclone/psyir/backend/sympy_writer.py
+++ b/src/psyclone/psyir/backend/sympy_writer.py
@@ -39,16 +39,21 @@
 '''
 
 import keyword
+from typing import Dict, List, Optional
 
 import sympy
 from sympy.parsing.sympy_parser import parse_expr
 
+from psyclone.core import (Signature, SingleVariableAccessInfo,
+                           VariablesAccessMap)
 from psyclone.psyir.backend.fortran import FortranWriter
 from psyclone.psyir.backend.visitor import VisitorError
 from psyclone.psyir.frontend.sympy_reader import SymPyReader
 from psyclone.psyir.nodes import (
-    DataNode, Range, Reference, IntrinsicCall, Call)
-from psyclone.psyir.symbols import (ArrayType, ScalarType, SymbolTable)
+    Call, DataNode, IntrinsicCall, Node, Range, Reference, StructureReference)
+from psyclone.psyir.symbols import (
+    ArrayType, DataSymbol, RoutineSymbol, ScalarType, Symbol,
+    SymbolError, SymbolTable, UnresolvedType)
 
 
 class SymPyWriter(FortranWriter):
@@ -109,10 +114,12 @@ class SymPyWriter(FortranWriter):
 
         # The writer will use special names in array expressions to indicate
         # the lower and upper bound (e.g. ``a(::)`` becomes
-        # ``a(sympy_lower, sympy_upper, 1)``). The symbol table will be used
-        # to resolve a potential name clash with a user variable.
+        # ``a(sympy_lower, sympy_upper, 1)`` while ``a`` becomes
+        # ``a(sympy_no_bounds, sympy_no_bounds, 1)``). The symbol table will
+        # be used to resolve a potential name clash with a user variable.
         self._lower_bound = "sympy_lower"
         self._upper_bound = "sympy_upper"
+        self._no_bounds = "sympy_no_bounds"
 
         if not SymPyWriter._RESERVED_NAMES:
             # Get the list of all reserved Python words from the
@@ -197,7 +204,11 @@ class SymPyWriter(FortranWriter):
                                   "never be called.")
 
     # -------------------------------------------------------------------------
-    def _create_sympy_array_function(self, name, sig=None, num_dims=None):
+    def _create_sympy_array_function(self,
+                                     name: str,
+                                     sig: Optional[Signature] = None,
+                                     num_dims: Optional[List[int]] = None,
+                                     is_call: Optional[bool] = False):
         '''Creates a Function class with the given name to be used for SymPy
         parsing. This Function overwrites the conversion to string, and will
         replace the triplicated array indices back to the normal Fortran
@@ -210,29 +221,29 @@ class SymPyWriter(FortranWriter):
         :param sig: the signature of the variable, which is required
             to convert user defined types back properly. Only defined for
             user-defined types.
-        :type sig: Optional[:py:class:`psyclone.core.Signature`]
         :param num_dims: the number of dimensions for each component of a
             user defined type.
-        :type num_dims: Optional[List[int]]
+        :param is_call: whether the Function represents an actual Call to
+            a function. If so, it is permitted to have no arguments.
 
         :returns: a SymPy function, which has a special ``_sympystr`` function
             defined as attribute to print user-defined types..
         :rtype: :py:class:`sympy.Function`
-        '''
 
-        # Now a new Fortran array is used. Create a new function
-        # instance, and overwrite how this function is converted back
-        # into a string by defining the ``_sympystr`` attribute,
-        # which points to a function that controls how this object
-        # is converted into a string. Use the ``print_fortran_array``
-        # function from the SymPyReader for this. Note that we cannot
-        # create a derived class based on ``Function`` and define
-        # this function there: SymPy tests internally if the type is a
-        # Function (not if it is an instance), therefore, SymPy's
-        # behaviour would change if we used a derived class:
+        '''
+        # Create a new sympy function instance, and overwrite how this
+        # function is converted back into a string by defining the
+        # ``_sympystr`` attribute, which points to a function that
+        # controls how this object is converted into a string. Use the
+        # ``print_fortran_array`` function from the SymPyReader for
+        # this. Note that we cannot create a derived class based on
+        # ``Function`` and define this function there: SymPy tests
+        # internally if the type is a Function (not if it is an
+        # instance), therefore, SymPy's behaviour would change if we
+        # used a derived class:
         # https://docs.sympy.org/latest/modules/functions/index.html:
-        # "It [Function class] also serves as a constructor for undefined
-        # function classes."
+        # "It [Function class] also serves as a constructor for
+        # undefined function classes."
         new_func = sympy.Function(name)
         # pylint: disable=protected-access
         new_func._sympystr = SymPyReader.print_fortran_array
@@ -242,13 +253,80 @@ class SymPyWriter(FortranWriter):
         # the indices back to the user defined types.
         new_func._sig = sig
         new_func._num_dims = num_dims
+        # We also have to store whether or not this is a function call in
+        # case it has no arguments.
+        new_func._is_call = is_call
         # pylint: enable=protected-access
         return new_func
 
+    @staticmethod
+    def _ndims_for_struct_access(sig: Signature,
+                                 sva: SingleVariableAccessInfo) -> List[int]:
+        '''
+        The same Signature can be accessed with different numbers of indices,
+        e.g. a%b, a%b(1) and a(1)%b. This routine examines all accesses and
+        returns a list of the maximum number of dimensions that is used for
+        each component of the Signature.
+
+        :param sig: the Signature of the structure access.
+        :param sva: details on all accesses for the Signature.
+
+        :returns: the number of dimensions of each component of the structure
+                  access represented by the supplied Signature.
+        '''
+        # This list will hold lists of the number of indices on
+        # each component, one list per access, e.g.:
+        #  a%b => [0,0] and  a(2)%b => [1,0]
+        num_dims_for_access = []
+        for access in sva.all_accesses:
+            indices = access.component_indices
+            # Create the list of number of indices on each component for
+            # this access.
+            num_dims = []
+            for i in range(len(sig)):
+                num_dims.append(len(indices[i]))
+            num_dims_for_access.append(num_dims)
+        # For each component, find the maximum number of dimensions
+        # seen in any access.
+        max_dims = []
+        for i in range(len(sig)):
+            max_dims.append(max(dims[i] for dims in num_dims_for_access))
+        return max_dims
+
+    @staticmethod
+    def _specialise_array_symbol(sym: Symbol, sva: SingleVariableAccessInfo):
+        '''
+        If we can be confident that the supplied Symbol should be of ArrayType
+        due to the way it is accessed then we specialise it in place.
+
+        :param sym: the Symbol to specialise.
+        :param sva: information on the ways in which the Symbol is accessed.
+
+        '''
+        if all(acs.is_array() for acs in sva.all_accesses):
+            return
+        if not sym or isinstance(sym, (DataSymbol, RoutineSymbol)):
+            return
+        # Find an access that has indices.
+        for acs in sva.all_accesses:
+            if not acs.is_array():
+                continue
+            ndims = None
+            for indices in acs.component_indices:
+                if indices:
+                    ndims = len(indices)
+            if ndims is not None:
+                sym.specialise(
+                    DataSymbol,
+                    datatype=ArrayType(UnresolvedType(),
+                                       [ArrayType.Extent.DEFERRED]*ndims))
+            return
+
     # -------------------------------------------------------------------------
-    def _create_type_map(self, list_of_expressions, identical_variables=None,
-                         all_variables_positive=None):
-        '''This function creates a dictionary mapping each Reference in any
+    def _create_type_map(self, list_of_expressions: List[Node],
+                         identical_variables: Optional[Dict[str, str]] = None,
+                         all_variables_positive: Optional[bool] = None):
+        '''This function creates a dictionary mapping each access in any
         of the expressions to either a SymPy Function (if the reference
         is an array reference) or a Symbol (if the reference is not an
         array reference). It defines a new SymPy function for each array,
@@ -256,6 +334,15 @@ class SymPyWriter(FortranWriter):
         converts array indices back by combining each three arguments into
         one expression (i. e. ``a(1,9,2)`` would become ``a(1:9:2)``).
 
+        An access like ``a(i)%b(j)`` is mapped to a function named ``a_b``
+        with arguments ``(i,i,1,j,j,1)`` (also handling name clashes in case
+        that the user code already contains a symbol ``a_b``). The SymPy
+        function created for this new symbol will store the original signature
+        and the number of indices for each member (so in the example above
+        that would be ``Signature("a%b")`` and ``(1,1)``. This information
+        is sufficient to convert the SymPy symbol back to the correct Fortran
+        representation.
+
         A new symbol table is created any time this function is called, so
         it is important to provide all expressions at once for the symbol
         table to avoid name clashes in any expression.
@@ -282,18 +369,15 @@ class SymPyWriter(FortranWriter):
         :param list_of_expressions: the list of expressions from which all
             references are taken and added to a symbol table to avoid
             renaming any symbols (so that only member names will be renamed).
-        :type list_of_expressions: List[:py:class:`psyclone.psyir.nodes.Node`]
         :param identical_variables: which variable names are known to represent
             identical quantities.
-        :type identical_variables: Optional[dict[str, str]]
-        :param Optional[bool] all_variables_positive: whether or not (the
-            default) to assume that all variables are positive definite
-            quantities.
+        :param all_variables_positive: whether or not (the default) to assume
+            that all variables are positive definite quantities.
 
         '''
-        # Create a new symbol table, so previous symbol will not affect this
+        # Create a new symbol table, so previous symbols will not affect this
         # new conversion (i.e. this avoids name clashes with a previous
-        # conversion). First add all reserved names so that these names will
+        # conversion). First, add all reserved names so that these names will
         # automatically be renamed. The symbol table is used later to also
         # create guaranteed unique names for lower and upper bounds.
         self._symbol_table = SymbolTable()
@@ -305,49 +389,67 @@ class SymPyWriter(FortranWriter):
         if all_variables_positive:
             assumptions["positive"] = True
 
-        # Find each reference in each of the expression, and declare this name
+        # Find every symbol in each of the expressions, and declare this name
         # as either a SymPy Symbol (scalar reference), or a SymPy Function
-        # (an array).
+        # (either an array or a function call).
+        vam = VariablesAccessMap()
         for expr in list_of_expressions:
-            # TODO #2542. References should be iterated with the
-            # reference_acess method when its issues are fixed.
-            for ref in expr.walk(Reference):
-                if (isinstance(ref.parent, Call) and
-                        ref.parent.children[0] is ref):
-                    continue
-
-                name = ref.name
-                # The reserved Python keywords do not have tags, so they
-                # will not be found.
-                if name in self._symbol_table.tags_dict:
-                    continue
-
-                # Any symbol from the list of expressions to be handled
-                # will be created with a tag, so if the same symbol is
-                # used more than once, the previous test will prevent
-                # calling new_symbol again. If the name is a Python
-                # reserved symbol, a new unique name will be created by
-                # the symbol table.
-                unique_sym = self._symbol_table.new_symbol(name, tag=name)
-                # Test if an array or an array expression is used:
-                if not ref.is_array:
-                    self._sympy_type_map[unique_sym.name] = sympy.Symbol(
-                        name, **assumptions)
-                    continue
+            vam.update(expr.reference_accesses())
+
+        for sig in vam.all_signatures:
+            sva: SingleVariableAccessInfo = vam[sig]
 
-                # A Fortran array is used which has not been seen before.
-                # Declare a new SymPy function for it. This SymPy function
-                # will convert array expressions back into the original
-                # Fortran code.
-                self._sympy_type_map[unique_sym.name] = \
-                    self._create_sympy_array_function(name)
+            flat_name = "_".join(name for name in sig)
+            unique_sym = self._symbol_table.find_or_create_tag(
+                str(sig), root_name=flat_name)
+
+            try:
+                # Depending on the situation, we won't always
+                # have a scope, hence the try...except.
+                orig_sym = sva.all_accesses[0].node.scope.symbol_table.\
+                    lookup(sig.var_name)
+            except SymbolError:
+                if isinstance(sva.all_accesses[0].node, Reference):
+                    orig_sym = sva.all_accesses[0].node.symbol
+                else:
+                    orig_sym = None
+
+            is_fn_call = isinstance(orig_sym, RoutineSymbol)
+
+            if (sva.is_array() or
+                    (orig_sym and (orig_sym.is_array or is_fn_call))):
+                # A Fortran array or function call. Declare a new SymPy
+                # function for it. This SymPy function will convert array
+                # expressions back into the original Fortran code.
+                if sig.is_structure:
+                    max_dims = self._ndims_for_struct_access(sig, sva)
+                    # We store the signature and the list of dimensions
+                    # associated with each component as part of the Sympy fn.
+                    self._sympy_type_map[unique_sym.name] = \
+                        self._create_sympy_array_function(
+                            unique_sym.name, sig, max_dims, is_call=is_fn_call)
+                else:
+                    # Not a structure access.
+                    self._sympy_type_map[unique_sym.name] = \
+                        self._create_sympy_array_function(sig.var_name,
+                                                          is_call=is_fn_call)
+                    # To avoid confusion in sympy_reader, we specialise any
+                    # Symbol that we are now confident is an array.
+                    self._specialise_array_symbol(orig_sym, sva)
+            else:
+                # A scalar access.
+                if sig.is_structure:
+                    self._sympy_type_map[unique_sym.name] = sympy.Symbol(
+                        sig.to_language(), **assumptions)
+                else:
+                    self._sympy_type_map[unique_sym.name] = sympy.Symbol(
+                        sig.var_name, **assumptions)
 
         if not identical_variables:
             identical_variables = {}
         # For all variables that are the same, set the symbols to be
         # identical. This means if e.g. identical_variables={'i': 'j'},
         # the expression i-j becomes j-j = 0
-
         for var1, var2 in identical_variables.items():
             if var1 in self._sympy_type_map and var2 in self._sympy_type_map:
                 self._sympy_type_map[var1] = self._sympy_type_map[var2]
@@ -360,25 +462,34 @@ class SymPyWriter(FortranWriter):
         self._upper_bound = \
             self._symbol_table.new_symbol("sympy_upper",
                                           tag="sympy!upper_bound").name
+        # This one is used when the array symbol is accessed (in the original
+        # code) without any indexing at all.
+        self._no_bounds = self._symbol_table.new_symbol(
+            "sympy_no_bounds", tag="sympy!no_bounds").name
 
     # -------------------------------------------------------------------------
     @property
-    def lower_bound(self):
+    def lower_bound(self) -> str:
         ''':returns: the name to be used for an unspecified lower bound.
-        :rtype: str
-
         '''
         return self._lower_bound
 
     # -------------------------------------------------------------------------
     @property
-    def upper_bound(self):
+    def upper_bound(self) -> str:
         ''':returns: the name to be used for an unspecified upper bound.
-        :rtype: str
-
         '''
         return self._upper_bound
 
+    # -------------------------------------------------------------------------
+    @property
+    def no_bounds(self) -> str:
+        '''
+        :returns: the name to be used when no bounds are present on an
+                  array access.
+        '''
+        return self._no_bounds
+
     # -------------------------------------------------------------------------
     @property
     def type_map(self):
@@ -536,58 +647,37 @@ class SymPyWriter(FortranWriter):
         return self.arrayofstructuresreference_node(node)
 
     # -------------------------------------------------------------------------
-    def arrayofstructuresreference_node(self, node):
-        '''This handles ArrayOfStructureReferences (and also simple
-        StructureReferences). An access like ``a(i)%b(j)`` is converted to
-        the string ``a_b(i,i,1,j,j,1)`` (also handling name clashes in case
-        that the user code already contains a symbol ``a_b``). The SymPy
-        function created for this new symbol will store the original signature
-        and the number of indices for each member (so in the example above
-        that would be ``Signature("a%b")`` and ``(1,1)``. This information
-        is sufficient to convert the SymPy symbol back to the correct Fortran
-        representation
+    def arrayofstructuresreference_node(self, node: StructureReference) -> str:
+        '''
+        This handles ArrayOfStructureReferences (and also simple
+        StructureReferences).
 
         :param node: a StructureReference PSyIR node.
-        :type node: :py:class:`psyclone.psyir.nodes.StructureReference`
 
-        :returns: the code as string.
-        :rtype: str
+        :returns: text representation of the code.
 
         '''
         sig, indices = node.get_signature_and_indices()
 
-        out = []
-        num_dims = []
         all_dims = []
-        is_array = False
         for i, name in enumerate(sig):
-            num_dims.append(len(indices[i]))
-            for index in indices[i]:
-                all_dims.append(index)
-                is_array = True
-            out.append(name)
-        flat_name = "_".join(out)
-
-        # Find (or create) a unique variable name:
-        try:
-            unique_name = self._symbol_table.lookup_with_tag(str(sig)).name
-        except KeyError:
-            unique_name = self._symbol_table.new_symbol(flat_name,
-                                                        tag=str(sig)).name
-        if is_array:
+            if indices[i]:
+                for index in indices[i]:
+                    all_dims.append(index)
+            else:
+                sym = self._symbol_table.lookup_with_tag("sympy!no_bounds")
+                all_dims.append(Reference(sym))
+
+        # Find the unique variable name created in _create_type_map()
+        unique_name = self._symbol_table.lookup_with_tag(str(sig)).name
+
+        if self.type_map[unique_name].is_Function:
+            # If the SymPy type is a Function then this means the original
+            # expression was either an array access or a function call.
             indices_str = self.gen_indices(all_dims)
-            # Create the corresponding SymPy function, which will store
-            # the signature and num_dims, so that the correct Fortran
-            # representation can be recreated later.
-            self._sympy_type_map[unique_name] = \
-                self._create_sympy_array_function(unique_name, sig, num_dims)
             return f"{unique_name}({','.join(indices_str)})"
 
-        # Not an array access. We use the unique name  for the string,
-        # but the required symbol is mapped to the original name, which means
-        # if the SymPy expression is converted to a string (in order to be
-        # parsed), it will use the original structure reference syntax:
-        self._sympy_type_map[unique_name] = sympy.Symbol(sig.to_language())
+        # Not an array access. We use the unique name for the string.
         return unique_name
 
     # -------------------------------------------------------------------------
@@ -621,6 +711,30 @@ class SymPyWriter(FortranWriter):
         # information can be ignored.
         return node.value
 
+    def call_node(self, node: Call) -> str:
+        '''
+        Handler for Call nodes in the PSyIR tree. Ensures that we use the
+        name created in the type map and encodes the arguments to the call
+        in the same form as is used for array accesses (since both are mapped
+        to a Sympy function).
+
+        Note that a :class:`Call` in the PSyIR can represent a call to either a
+        subroutine or a function. The latter can appear within expressions and
+        either can have zero arguments.
+
+        :param node: a PSyIR Call node.
+
+        :returns: the SymPy representation for the Call.
+
+        '''
+        target_name = node.routine.symbol.name
+        # Find the unique variable name created in _create_type_map()
+        unique_name = self._symbol_table.lookup_with_tag(target_name).name
+        # Encode the arguments to the call (if any) in the same way as we do
+        # for indices on an array access.
+        indices_str = self.gen_indices(node.arguments)
+        return f"{unique_name}({','.join(indices_str)})"
+
     def intrinsiccall_node(self, node):
         ''' This method is called when an IntrinsicCall instance is found in
         the PSyIR tree. The Sympy backend will use the exact sympy name for
@@ -665,7 +779,7 @@ class SymPyWriter(FortranWriter):
         PSyIR tree. It handles the case that this normal reference might
         be an array expression, which in the SymPy writer needs to have
         indices added explicitly: it basically converts the array expression
-        ``a`` to ``a(sympy_lower, sympy_upper, 1)``.
+        ``a`` to ``a(sympy_no_bounds, sympy_no_bounds, 1)``.
 
         :param node: a Reference PSyIR node.
         :type node: :py:class:`psyclone.psyir.nodes.Reference`
@@ -683,16 +797,17 @@ class SymPyWriter(FortranWriter):
             # If the tag did not exist it means that this symbol has not
             # been re-named, and we can use it as is.
             name = node.name
+
         if not node.is_array:
             # This reference is not an array, just return the name
             return name
 
-        # Now this must be an array expression without parenthesis. Add
-        # the triple-array indices to represent `lower:upper:1` for each
+        # Now this must be an array expression without parentheses. For
+        # consistency, we still treat it as a Sympy function call and therefore
+        # add the triple array indices to represent `lower:upper:1` for each
         # dimension:
         shape = node.symbol.shape
-        result = [f"{self.lower_bound},"
-                  f"{self.upper_bound},1"]*len(shape)
+        result = [f"{self.no_bounds},{self.no_bounds},1"]*len(shape)
 
         return (f"{name}{self.array_parenthesis[0]}"
                 f"{','.join(result)}{self.array_parenthesis[1]}")
diff --git a/src/psyclone/psyir/frontend/fparser2.py b/src/psyclone/psyir/frontend/fparser2.py
index 465232898..0671b985c 100644
--- a/src/psyclone/psyir/frontend/fparser2.py
+++ b/src/psyclone/psyir/frontend/fparser2.py
@@ -4456,10 +4456,12 @@ class Fparser2Reader():
             # TODO 2884: We should be able to handle this imported symbol
             # better. If we can, we need to handle a case where is_elemental
             # can be None.
-            if isinstance(ref.symbol.interface, ImportInterface):
+            if isinstance(ref.symbol.interface, (ImportInterface,
+                                                 UnresolvedInterface)):
                 raise NotImplementedError(
-                        "PSyclone doesn't yet support reference to imported "
-                        "symbols inside WHERE clauses.")
+                        f"PSyclone doesn't yet support references to imported/"
+                        f"unresolved symbols inside WHERE clauses: "
+                        f"'{ref.symbol.name}' is unresolved.")
             if (isinstance(ref.symbol, DataSymbol) and
                     elemental_ancestor):
                 try:
@@ -4477,7 +4479,7 @@ class Fparser2Reader():
             raise NotImplementedError(
                 f"Only WHERE constructs using explicit array notation "
                 f"including ranges (e.g. 'my_array(1,:)') are supported but "
-                f"found '{logical_expr}'")
+                f"found '{logical_expr[0]}'")
 
         array_ref = first_array.ancestor(Reference, include_self=True)
         if not isinstance(array_ref.datatype, ArrayType):
diff --git a/src/psyclone/psyir/frontend/sympy_reader.py b/src/psyclone/psyir/frontend/sympy_reader.py
index 912636e71..fc3cd8520 100644
--- a/src/psyclone/psyir/frontend/sympy_reader.py
+++ b/src/psyclone/psyir/frontend/sympy_reader.py
@@ -91,10 +91,12 @@ class SymPyReader():
     # correct values in the constructor.
     _lower_bound = "sympy_lower"
     _upper_bound = "sympy_upper"
+    _no_bounds = "no_sympy_bounds"
 
     def __init__(self, sympy_writer):
         SymPyReader._lower_bound = sympy_writer.lower_bound
         SymPyReader._upper_bound = sympy_writer.upper_bound
+        SymPyReader._no_bounds = sympy_writer.no_bounds
 
     # -------------------------------------------------------------------------
     def psyir_from_expression(self, sympy_expr, symbol_table):
@@ -125,7 +127,7 @@ class SymPyReader():
 
     # -------------------------------------------------------------------------
     # pylint: disable=no-self-argument, too-many-branches
-    def print_fortran_array(function, printer):
+    def print_fortran_array(function, printer) -> str:
         '''A custom print function to convert a modified Fortran array access
         back to standard Fortran. This function is set as ``_sympystr_`` method
         of the SymPy functions created in the SymPyWriter (see
@@ -134,9 +136,11 @@ class SymPyReader():
         convert being the first argument! This function converts
         the three values that each index is converted to back into the Fortran
         array notation. It uses the class variables
-        ``SymPyReader._lower_bound`` and ``SymPyReader._upper_bound`` as the
-        names that were used when the SymPy expressions were created in order
-        to convert array expressions correctly back.
+        ``SymPyReader._no_bounds``, ``SymPyReader._lower_bound`` and
+        ``SymPyReader._upper_bound`` as the names that were used when the
+        SymPy expressions were originally created in order
+        to correctly convert back array expressions (including the case where
+        the original array reference did not have parentheses).
 
         :param function: this function is called from a SymPy Function class,
             therefore the first argument is a SymPy Function instance (and NOT
@@ -146,7 +150,6 @@ class SymPyReader():
         :type printer: :py:class:`sympy.printing.str.StrPrinter`
 
         :returns: the string representation of this array access.
-        :rtype: str
 
         '''
         # pylint: disable=protected-access, no-member
@@ -154,11 +157,16 @@ class SymPyReader():
         name = function.__class__.__name__
         lower_b = SymPyReader._lower_bound
         upper_b = SymPyReader._upper_bound
+        no_bounds = SymPyReader._no_bounds
 
         # Analyse each triple of parameters, and add the corresponding
         # converted index (or array expression) to new_args:
         new_args = []
         for i in range(0, len(args), 3):
+            if args[i] == no_bounds and args[i+1] == no_bounds:
+                # This reference didn't have parentheses originally so we
+                # don't add any range information.
+                continue
             if args[i] == args[i+1] and args[i+2] == "1":
                 # a(i,i,1) --> a(i)
                 new_args.append(args[i])
@@ -189,8 +197,14 @@ class SymPyReader():
                     new_args.append(f"{args[i]}:{args[i+1]}:"
                                     f"{args[i+2]}")
 
-        if function._sig is None:
-            # It's not a user defined type, just create the array access:
+        if not new_args and not function._is_call:
+            # This reference is not a function call and didn't have parentheses
+            # originally so we reproduce that input.
+            return name
+
+        if function._sig is None or function._is_call:
+            # It's not a user defined type, just create the array access
+            # or function call:
             return f"{name}({','.join(new_args)})"
 
         # It is a user defined type. Re-assemble the original call by
diff --git a/src/psyclone/psyir/transformations/inline_trans.py b/src/psyclone/psyir/transformations/inline_trans.py
index 3ac3dc977..9c19cd397 100644
--- a/src/psyclone/psyir/transformations/inline_trans.py
+++ b/src/psyclone/psyir/transformations/inline_trans.py
@@ -187,7 +187,7 @@ class InlineTrans(Transformation):
         for child in node.arguments:
             try:
                 # TODO #1858, this won't yet work for arrays inside structures.
-                ref2arraytrans.apply(child)
+                ref2arraytrans.apply(child, allow_call_arguments=True)
             except (TransformationError, ValueError):
                 pass
 
diff --git a/src/psyclone/psyir/transformations/reference2arrayrange_trans.py b/src/psyclone/psyir/transformations/reference2arrayrange_trans.py
index 9786671ea..f544d244c 100644
--- a/src/psyclone/psyir/transformations/reference2arrayrange_trans.py
+++ b/src/psyclone/psyir/transformations/reference2arrayrange_trans.py
@@ -43,13 +43,15 @@
 '''
 from psyclone.errors import LazyString
 from psyclone.psyGen import Transformation
-from psyclone.psyir.nodes import (Range, Reference, ArrayReference, Literal,
-                                  IntrinsicCall, Assignment)
-from psyclone.psyir.symbols import INTEGER_TYPE, ArrayType
-from psyclone.psyir.transformations.transformation_error \
-    import TransformationError
+from psyclone.psyir.nodes import (ArrayReference, Assignment, Call,
+                                  IntrinsicCall, Literal, Range, Reference)
+from psyclone.psyir.symbols import INTEGER_TYPE, ArrayType, Symbol
+from psyclone.psyir.transformations.transformation_error import (
+    TransformationError)
+from psyclone.utils import transformation_documentation_wrapper
 
 
+@transformation_documentation_wrapper
 class Reference2ArrayRangeTrans(Transformation):
     '''Provides a transformation from PSyIR Array Notation (a reference to
     an Array) to a PSyIR Range. For example:
@@ -124,22 +126,28 @@ class Reference2ArrayRangeTrans(Transformation):
         step = Literal("1", INTEGER_TYPE)
         return (lower_bound, upper_bound, step)
 
-    def validate(self, node, options=None):
+    def validate(self, node, **kwargs):
         '''Check that the node is a Reference node and that the symbol it
         references is an array.
 
         :param node: a Reference node.
         :type node: :py:class:`psyclone.psyir.nodes.Reference`
-        :param options: a dict with options for transformations.
-        :type options: Optional[Dict[str, Any]]
+        :param allow_call_arguments: by default, any references that may be
+            arguments to non-elemental routines are not transformed. However,
+            this transformation is sometimes used in other transformations
+            where this restriction does not apply.
 
         :raises TransformationError: if the node is not a Reference
             node or the Reference node not does not reference an array
             symbol.
-        :raises TransformationError: if the Reference node is within an
-            inquiry or DEALLOCATE intrinsic.
+        :raises TransformationError: if the Reference node is (or may be)
+            passed as an argument to a call that is not elemental and
+            `allow_call_arguments` is False.
 
         '''
+        self.validate_options(**kwargs)
+        allow_call_arguments = self.get_option("allow_call_arguments",
+                                               **kwargs)
         # TODO issue #1858. Add support for structures containing arrays.
         # pylint: disable=unidiomatic-typecheck
         if not type(node) is Reference:
@@ -150,17 +158,22 @@ class Reference2ArrayRangeTrans(Transformation):
             raise TransformationError(
                 f"The supplied node should be a Reference to a symbol "
                 f"that is an array, but '{node.symbol.name}' is not.")
-        if isinstance(node.parent, IntrinsicCall) and node.parent.is_inquiry:
-            raise TransformationError(
-                f"References to arrays passed as arguments to intrinsic "
-                f"enquiry routine '{node.parent.routine.name}' should not be "
-                f"transformed.")
-        if (isinstance(node.parent, IntrinsicCall) and
-                node.parent.routine.name in ["DEALLOCATE"]):
+        if not allow_call_arguments and (isinstance(node.parent, Call) and
+                                         not node.parent.is_elemental):
+            raise TransformationError(LazyString(
+                lambda: f"The supplied node is passed as an argument to a "
+                f"Call to a non-elemental routine ("
+                f"{node.parent.debug_string().strip()}) and should not be "
+                f"transformed."))
+        if (isinstance(node.parent, Reference) and (
+                type(node.parent.symbol) is Symbol
+                or not isinstance(node.parent.symbol.datatype, ArrayType))):
             raise TransformationError(LazyString(
-                lambda: f"References to arrays passed to "
-                f"'{node.parent.routine.name}' intrinsics should not be "
-                f"transformed, but found:\n {node.parent.debug_string()}"))
+                lambda: f"References to arrays that *may* be routine arguments"
+                f" should not be transformed but found:\n "
+                f"{node.parent.debug_string()} and {node.parent.symbol.name} "
+                f"is not known to be of ArrayType (and therefore may be a "
+                f"call)."))
         assignment = node.ancestor(Assignment)
         if assignment and assignment.is_pointer:
             raise TransformationError(
@@ -168,7 +181,7 @@ class Reference2ArrayRangeTrans(Transformation):
                 f" inside pointer assignments, but found '{node.name}' in"
                 f" {assignment.debug_string()}")
 
-    def apply(self, node, options=None):
+    def apply(self, node, allow_call_arguments: bool = False, **kwargs):
         '''Apply the Reference2ArrayRangeTrans transformation to the specified
         node. The node must be a Reference to an array. The Reference
         is replaced by an ArrayReference with appropriate explicit
@@ -176,11 +189,13 @@ class Reference2ArrayRangeTrans(Transformation):
 
         :param node: a Reference node.
         :type node: :py:class:`psyclone.psyir.nodes.Reference`
-        :param options: a dict with options for transformations.
-        :type options: Optional[Dict[str, Any]]
+        :param allow_call_arguments: by default, any references that may be
+            arguments to non-elemental routines are not transformed. However,
+            this transformation is sometimes used in other transformations
+            where this restriction does not apply.
 
         '''
-        self.validate(node, options=None)
+        self.validate(node, allow_call_arguments=allow_call_arguments)
 
         symbol = node.symbol
         indices = []


Test patch:
diff --git a/src/psyclone/tests/core/symbolic_maths_test.py b/src/psyclone/tests/core/symbolic_maths_test.py
index 00d96a3c3..7b7ce82a1 100644
--- a/src/psyclone/tests/core/symbolic_maths_test.py
+++ b/src/psyclone/tests/core/symbolic_maths_test.py
@@ -42,6 +42,7 @@ from sympy import solvers, Symbol
 
 from psyclone.core.symbolic_maths import SymbolicMaths
 from psyclone.psyir.backend.sympy_writer import SymPyWriter
+from psyclone.psyir.nodes import Assignment
 
 
 def test_sym_maths_get():
@@ -122,7 +123,6 @@ def test_symbolic_math_equal(fortran_reader, expressions):
 
 
 @pytest.mark.parametrize("expressions", [("a%b", "a%b"),
-                                         ("c", "c(::,::)"),
                                          ("a%b(i)", "a%b(i)"),
                                          ("a%b(2*i)", "a%b(3*i-i)"),
                                          ("a%b(i-1)%c(j+1)",
@@ -475,7 +475,12 @@ def test_symbolic_math_use_range(fortran_reader, expressions):
     ("a*(b+c)", "a * b + a * c"),
     ("a*((b+c)/d)", "a * b / d + a * c / d"),
     ("a(i)*((b(i,j)+c(j))/d)",
-     "a(i) * b(i,j) / d + a(i) * c(j) / d")])
+     "a(i) * b(i,j) / d + a(i) * c(j) / d"),
+    # 'a' is unresolved so we don't know from the first occurrence whether or
+    # not it is a scalar.
+    ("a / a(i)", "a / a(i)"),
+    ("norm_u(idx+iw2) * u_e(idx + (LBOUND(u_e,dim=1)-iw2v), df2)",
+     "norm_u(idx + iw2) * u_e(idx - iw2v + LBOUND(u_e, 1),df2)")])
 def test_symbolic_maths_expand(fortran_reader, fortran_writer, expr, expected):
     '''Test the expand method works as expected.'''
     # A dummy program to easily create the PSyIR for the
@@ -492,6 +497,93 @@ def test_symbolic_maths_expand(fortran_reader, fortran_writer, expr, expected):
     assert result == expected
 
 
+def test_expand_with_intrinsic(fortran_reader, fortran_writer):
+    '''
+    Test that calling the `expand` method does not alter array accesses
+    that are passed as arguments to other routines - in this case the
+    LBOUND intrinsic.
+
+    '''
+    source = '''
+  subroutine apply_mixed_operator_code(ncell1, nlayers, ndf_w2, ndf_w2h, &
+                                  undf_w2v, undf_w2, lhs_w, norm_u, mu_cd)
+    use kinds_mod, only: i_def, r_solver
+    integer, intent(in) :: nlayers, ndf_w2, ndf_w2h, ncell1
+    integer, intent(in) :: undf_w2v, undf_w2
+    real(kind=r_solver), dimension(undf_w2v), intent(inout) :: lhs_w
+    real(kind=r_solver), dimension(undf_w2), intent(in) :: norm_u
+    real(kind=r_solver), dimension(ncell1,ndf_w2,ndf_w2), intent(in) :: mu_cd
+    real(kind=r_solver), dimension(0:nlayers - 1,ndf_w2) :: u_e
+    real(kind=r_solver), dimension(0:nlayers) :: t_col
+    integer(kind=i_def) :: idx_10, df, df2, ij, iwt, iw2, iw2h, iw2v
+    lhs_w(idx_10) = norm_u(idx_10 + (iw2 - iw2v)) * &
+        mu_cd(idx_10 + (ij - iw2v), ndf_w2h + df ,df2) * &
+        u_e(idx_10 + (LBOUND(u_e, dim=1) - iw2v), df2)
+  end subroutine apply_mixed_operator_code'''
+    psyir = fortran_reader.psyir_from_source(source)
+    sym_maths = SymbolicMaths.get()
+    rhs = psyir.walk(Assignment)[0].rhs
+    sym_maths.expand(rhs)
+    result = fortran_writer(psyir).lower()
+    # Check that the 'u_e' argument remains unchanged.
+    assert "lbound(u_e, 1),df2)" in result
+
+
+def test_symbolic_maths_expand_function(fortran_reader, fortran_writer):
+    '''Test the expand method works as expected when one of the
+    Symbols in the expression corresponds to an actual function.
+
+    '''
+    source = (
+        "module test\n"
+        "  implicit none\n"
+        "  contains\n"
+        "  subroutine use_a\n"
+        "    use some_mod, only: b, c, d\n"
+        "    integer :: i, j, x\n"
+        "    x = a(i)*((b(i,j)+c(j))/d)\n"
+        "  end subroutine\n"
+        "  function a(i)\n"
+        "     integer, intent(in) :: i\n"
+        "     integer :: a\n"
+        "     a = i\n"
+        "  end function\n"
+        "end module\n")
+    psyir = fortran_reader.psyir_from_source(source)
+    sym_maths = SymbolicMaths.get()
+    assigns = psyir.walk(Assignment)
+    sym_maths.expand(assigns[0].rhs)
+    result = fortran_writer(psyir).lower()
+    assert "a(i) * b(i,j) / d +" in result
+
+
+def test_symbolic_maths_expand_function_no_arg(fortran_reader, fortran_writer):
+    '''Test the expand method works as expected when one of the Symbols in the
+    expression corresponds to a function call with no arguments.
+
+    '''
+    source = (
+        "module test\n"
+        "  implicit none\n"
+        "  contains\n"
+        "  subroutine use_a\n"
+        "    use some_mod, only: b, c, d\n"
+        "    integer :: x, i, j\n"
+        "    x = a()*((b(i,j)+c(j))/d)\n"
+        "  end subroutine\n"
+        "  function a()\n"
+        "     integer :: a\n"
+        "     a = 10\n"
+        "  end function\n"
+        "end module\n")
+    psyir = fortran_reader.psyir_from_source(source)
+    sym_maths = SymbolicMaths.get()
+    assigns = psyir.walk(Assignment)
+    sym_maths.expand(assigns[0].rhs)
+    result = fortran_writer(psyir).lower()
+    assert "x = a() * b(i,j) / d + a() *" in result
+
+
 def test_symbolic_maths_array_and_array_index(fortran_reader):
     '''Test having an expression that uses a whole array and
     the same array with an index, e.g. : `a(i) + a`.
@@ -503,12 +595,10 @@ def test_symbolic_maths_array_and_array_index(fortran_reader):
           y = a
         end program test_prog'''
     psyir = fortran_reader.psyir_from_source(source)
+    assigns = psyir.walk(Assignment)
     sym_maths = SymbolicMaths.get()
-    assert not sym_maths.equal(psyir.children[0][0].rhs,
-                               psyir.children[0][1].rhs)
+    assert not sym_maths.equal(assigns[0].rhs, assigns[1].rhs)
 
-    assert sym_maths.equal(psyir.children[0][0].rhs,
-                           psyir.children[0][0].rhs)
+    assert sym_maths.equal(assigns[0].rhs, assigns[0].rhs)
 
-    assert sym_maths.equal(psyir.children[0][1].rhs,
-                           psyir.children[0][1].rhs)
+    assert sym_maths.equal(assigns[1].rhs, assigns[1].rhs)
diff --git a/src/psyclone/tests/psyir/backend/sympy_writer_test.py b/src/psyclone/tests/psyir/backend/sympy_writer_test.py
index 99572651f..60a574aa8 100644
--- a/src/psyclone/tests/psyir/backend/sympy_writer_test.py
+++ b/src/psyclone/tests/psyir/backend/sympy_writer_test.py
@@ -165,7 +165,8 @@ def test_sym_writer_real_constants(fortran_reader, expressions):
                                          ("MIN(1,2)", "Min(1, 2)"),
                                          ("MOD(1,2)", "Mod(1, 2)"),
                                          ("EXP(1)", "exp(1)"),
-                                         ("LBOUND(1,2)", "LBOUND(1, 2)")
+                                         ("LBOUND(1,2)", "LBOUND(1, 2)"),
+                                         ("LBOUND(1,dim=2)", "LBOUND(1, 2)")
                                          ])
 def test_sym_writer_functions(fortran_reader, expressions):
     '''Test that certain standard functions are recognised and
@@ -186,26 +187,29 @@ def test_sym_writer_functions(fortran_reader, expressions):
     assert SymPyWriter()._to_str(function) == expressions[1]
 
 
-@pytest.mark.parametrize("expr, sym_map", [("i", {'i': Symbol('i')}),
-                                           ("f(1)", {'f': Function('f')}),
-                                           ("f(:)", {'f': Function('f')}),
-                                           ("a%b", {'a': Symbol('a'),
-                                                    'a_b': Symbol('a_b')}),
-                                           ("a%b(1)", {'a': Symbol('a'),
-                                                       'a_b': Function('a_b')})
-                                           ])
+@pytest.mark.parametrize("expr, sym_map",
+                         [("i", {'i': Symbol('i')}),
+                          ("f(1)", {'f': Function('f')}),
+                          ("f(:)", {'f': Function('f')}),
+                          ("a%b", {'a_b': Symbol('a_b')}),
+                          ("a%b(1)", {'a_b': Function('a_b')}),
+                          ("c + c(1)", {'c': Function('c')}),
+                          ("a%b + a%b(1)", {'a_b': Function('a_b')}),
+                          # iptr will be of UnknownFortranType
+                          ("LBOUND(iptr)", {'iptr': Function('iptr')})
+                          ])
 def test_sympy_writer_create_type_map(expr, sym_map, fortran_reader):
     '''Tests that the static create_type_map creates a dictionary
-    with correctly declared references (and not any member names,
-    which will be added later).
-    '''
+    with correctly declared references.
 
+    '''
     # A dummy program to easily create the PSyIR for the
     # expressions we need. We just take the RHS of the assignments
     source = f'''program test_prog
                 use my_mod
                 type(my_type) :: a, b(10)
                 integer :: i, f(10)
+                integer, pointer :: iptr(5)
                 x = {expr}
                 end program test_prog '''
 
@@ -216,21 +220,22 @@ def test_sympy_writer_create_type_map(expr, sym_map, fortran_reader):
     assert sympy_writer._sympy_type_map.keys() == sym_map.keys()
 
 
-@pytest.mark.parametrize("expressions", [("a%x", "a_x"),
-                                         ("b(i)%x", "b_x(i,i,1)"),
-                                         ("a%x(i)", "a_x(i,i,1)"),
-                                         ("b(j)%x(i)", "b_x(j,j,1,i,i,1)"),
-                                         ("b(i)%c(b_c)",
-                                          "b_c_1(i,i,1,b_c,b_c,1)"),
-                                         ("a_c + a%c(i)",
-                                          "a_c + a_c_1(i,i,1)"),
-                                         ("b(b_c)%c(i)",
-                                          "b_c_1(b_c,b_c,1,i,i,1)"),
-                                         ("b(b_c)%c(i)",
-                                          "b_c_1(b_c,b_c,1,i,i,1)"),
-                                         ("a_b_c + a_b_c_1 + a%b%c",
-                                          "a_b_c + a_b_c_1 + a_b_c_2"),
-                                         ])
+@pytest.mark.parametrize(
+    "expressions",
+    [("a%x", "a_x"),
+     ("b(i)%x", "b_x(i,i,1,sympy_no_bounds,sympy_no_bounds,1)"),
+     ("a%x(i)", "a_x(sympy_no_bounds,sympy_no_bounds,1,i,i,1)"),
+     ("b(j)%x(i)", "b_x(j,j,1,i,i,1)"),
+     ("b(i)%c(b_c)", "b_c(i,i,1,b_c_1,b_c_1,1)"),
+     ("a_c + a%c(i)", "a_c_1 + a_c(sympy_no_bounds,sympy_no_bounds,1,i,i,1)"),
+     ("b(b_c)%c(i)", "b_c(b_c_1,b_c_1,1,i,i,1)"),
+     ("b(b_c)%c(i)", "b_c(b_c_1,b_c_1,1,i,i,1)"),
+     ("a_b_c + a_b_c_1 + a%b%c",
+      "a_b_c_1 + a_b_c_1_1 + a_b_c"),
+     ("b(1)%c + b(j)%c(1) + b(i)%c(2)",
+      "b_c(1,1,1,sympy_no_bounds,sympy_no_bounds,1) + b_c(j,j,1,1,1,1) + "
+      "b_c(i,i,1,2,2,1)")
+     ])
 def test_sym_writer_rename_members(fortran_reader, expressions):
     '''Test that members are converted and get a unique name that
     does not clash with any other variable used in the expression.
@@ -250,26 +255,25 @@ def test_sym_writer_rename_members(fortran_reader, expressions):
     assert SymPyWriter()._to_str(expr) == expressions[1]
 
 
-@pytest.mark.parametrize("expr, sym_map", [("a%x", {"a": Symbol("a"),
-                                                    "a_x": Symbol("a%x")}),
-                                           ("a%x(i)", {"a": Symbol("a"),
-                                                       "a_x": Function("a_x"),
-                                                       "i": Symbol("i")}),
-                                           ("b(i)%x(i)",
-                                            {"b": Function("b"),
-                                             "b_x": Function("b_x"),
-                                             "i": Symbol("i")}),
-                                           ("b(b_c)%c(i)",
-                                            {"b": Function("b"),
-                                             "b_c": Symbol("b_c"),
-                                             "b_c_1": Function("b_c_1"),
-                                             "i": Symbol("i")}),
-                                           ])
-def test_sym_writer_symbol_types(fortran_reader, expr, sym_map):
+@pytest.mark.parametrize(
+    "expr, positive, sym_map",
+    [("a%x", False, {"a_x": Symbol("a%x")}),
+     ("a%x", True, {"a_x": Symbol("a%x", **{"positive": True})}),
+     ("a%x(i)", False, {"a_x": Function("a_x"), "i": Symbol("i")}),
+     ("a%x(i)", True, {"a_x": Function("a_x"),
+                       "i": Symbol("i", **{"positive": True})}),
+     ("b(i)%x(i)", False, {"b_x": Function("b_x"), "i": Symbol("i")}),
+     ("b(b_c)%c(i)", False, {"b_c_1": Symbol("b_c"), "b_c": Function("b_c"),
+                             "i": Symbol("i")}),
+     ("b(b_c)%c(i)", True, {"b_c_1": Symbol("b_c", **{"positive": True}),
+                            "b_c": Function("b_c"),
+                            "i": Symbol("i", **{"positive": True})}),
+     ])
+def test_sym_writer_symbol_types(fortran_reader, expr, positive, sym_map):
     '''Tests that arrays are detected as SymPy functions, and scalars
-    as SymPy symbols. The expressions parameter contains as first
-    element the expression to parse, and as second element the
-    expected mapping of names to SymPy functions or symbols.
+    as SymPy symbols. The 'expr' parameter contains the expression to parse,
+    'positive' is whether or not to flag symbols as positive definite and
+    'sym_map' is the expected mapping of names to SymPy functions or symbols.
 
     '''
     # A dummy program to easily create the PSyIR for the
@@ -284,23 +288,21 @@ def test_sym_writer_symbol_types(fortran_reader, expr, sym_map):
     psyir = fortran_reader.psyir_from_source(source)
     expr = psyir.children[0].children[0].rhs
     sympy_writer = SymPyWriter()
-    _ = sympy_writer(expr)
-    assert sympy_writer.type_map.items() == sym_map.items()
-
-
-@pytest.mark.parametrize("expr, sym_map", [("i", {'i': Symbol('i')}),
-                                           ("f(1)", {'f': Function('f')}),
-                                           ("a%b", {'a': Symbol('a'),
-                                                    'a_b': Symbol('a_b')}),
-                                           ("a%b(1)",
-                                            {'a': Symbol('a'),
-                                             'a_b': Function('a_b')})
-                                           ])
-def test_sympy_writer_get_symbol_and_map(expr, sym_map, fortran_reader):
-    '''Tests that `get_sympy_expressions_and_symbol_map` function
-    returns the right symbol map.
-    '''
-
+    _ = sympy_writer(expr, all_variables_positive=positive)
+    assert len(sympy_writer.type_map) == len(sym_map)
+    for key in sympy_writer.type_map.keys():
+        assert sympy_writer.type_map[key] == sym_map[key]
+
+
+@pytest.mark.parametrize("expr, sym_map",
+                         [("i", {'i': Symbol('i')}),
+                          ("f(1)", {'f': Function('f')}),
+                          ("a%b", {'a_b': Symbol('a_b')}),
+                          ("a%b(1)", {'a_b': Function('a_b')}),
+                          ("lbound(f)", {'f': Function('f')}),
+                          ("1_i_def", {'i_def': Symbol('i_def')})])
+def test_sympy_writer_type_map(expr, sym_map, fortran_reader):
+    '''Tests that `_sympy_type_map` is setup correctly.'''
     # A dummy program to easily create the PSyIR for the
     # expressions we need. We just take the RHS of the assignments
     source = f'''program test_prog
@@ -317,11 +319,14 @@ def test_sympy_writer_get_symbol_and_map(expr, sym_map, fortran_reader):
     # Ignore the converted expressions here, they are tested elsewhere
     _ = writer([expr])
     assert writer._sympy_type_map.keys() == sym_map.keys()
+    # Repeat when the supplied Node does not have a parent tree
+    expr.detach()
+    _ = writer([expr])
+    assert writer._sympy_type_map.keys() == sym_map.keys()
 
 
-def test_sym_writer_convert_to_sympy_expressions(fortran_reader):
-    '''Tests that convenience function `convert_to_sympy_expressions`
-    works as expected.
+def test_sym_writer_parse_expr(fortran_reader):
+    '''Tests that convenience function `parse_expr` works as expected.
 
     '''
     # A dummy program to easily create the PSyIR for the
@@ -340,9 +345,10 @@ def test_sym_writer_convert_to_sympy_expressions(fortran_reader):
     sympy_writer = SymPyWriter()
     sympy_list = sympy_writer([exp1, exp2])
 
-    expr = parse_expr("a_b_1 + a_c(1,1,1) + i", sympy_writer.type_map)
+    expr = parse_expr("a_b + a_c(sympy_no_bounds,sympy_no_bounds,1,1,1,1) + i",
+                      sympy_writer.type_map)
     assert sympy_list[0] == expr
-    assert sympy_list[1] == parse_expr("a_b + j", sympy_writer.type_map)
+    assert sympy_list[1] == parse_expr("a_b_1 + j", sympy_writer.type_map)
 
 
 def test_sym_writer_parse_errors(fortran_reader):
@@ -368,60 +374,48 @@ def test_sym_writer_parse_errors(fortran_reader):
             in str(err.value))
 
 
-@pytest.mark.parametrize("expressions", [("b(i)", "b(i,i,1)"),
-                                         ("b(:)",
-                                          "b(sympy_lower,sympy_upper,1)"),
-                                         ("b(::)",
-                                          "b(sympy_lower,sympy_upper,1)"),
-                                         ("b(5::)", "b(5,sympy_upper,1)"),
-                                         ("b(:5:)", "b(sympy_lower,5,1)"),
-                                         ("b(::5)",
-                                          "b(sympy_lower,sympy_upper,5)"),
-                                         ("b(i::)", "b(i,sympy_upper,1)"),
-                                         ("b(:i:)", "b(sympy_lower,i,1)"),
-                                         ("b(::i)",
-                                          "b(sympy_lower,sympy_upper,i)"),
-                                         ("b(i:5:)", "b(i,5,1)"),
-                                         ("b(i:j:)", "b(i,j,1)"),
-                                         ("b(i::j)", "b(i,sympy_upper,j)"),
-                                         ("b(:i:j)", "b(sympy_lower,i,j)"),
-                                         ("b(i:j:k)", "b(i,j,k)"),
-                                         ("b", "b(sympy_lower,sympy_upper,1)"),
-                                         ("c(i,j)", "c(i,i,1,j,j,1)"),
-                                         ("c(::,::)",
-                                          "c(sympy_lower,sympy_upper,1,"
-                                          "sympy_lower,sympy_upper,1)"),
-                                         ("c", "c(sympy_lower,sympy_upper,1,"
-                                               "sympy_lower,sympy_upper,1)"),
-                                         ("b(i)%x", "b_x(i,i,1)"),
-                                         ("b(i)%x(j)", "b_x(i,i,1,j,j,1)"),
-                                         ("c(i,j)%x", "c_x(i,i,1,j,j,1)"),
-                                         ("c(i,j)%x(j)",
-                                          "c_x(i,i,1,j,j,1,j,j,1)"),
-                                         ("c(i,j)%d%e",
-                                          "c_d_e(i,i,1,j,j,1)"),
-                                         ("c(i,j)%d%f(i)",
-                                          "c_d_f(i,i,1,j,j,1,i,i,1)"),
-                                         ("c(i::k,j)%d%f(i:j:k)",
-                                          "c_d_f(i,sympy_upper,k,j,j,1,"
-                                          "i,j,k)"),
-                                         # Check name clashes, if a user
-                                         # variable is the same as the names
-                                         # for upper/lower bound
-                                         ("sympy_upper(:)",
-                                          "sympy_upper(sympy_lower,"
-                                          "sympy_upper_1,1)"),
-                                         ("sympy_lower(:)",
-                                          "sympy_lower(sympy_lower_1,"
-                                          "sympy_upper,1)"),
-                                         # The +sympy_upper at the end is
-                                         # an array expression, so it gets
-                                         # indices added!
-                                         ("sympy_lower(:)+sympy_upper",
-                                          "sympy_lower(sympy_lower_1,"
-                                          "sympy_upper_1,1) + sympy_upper"
-                                          "(sympy_lower_1,sympy_upper_1,1)"),
-                                         ])
+@pytest.mark.parametrize(
+    "expressions",
+    [("b(i)", "b(i,i,1)"),
+     ("b(:)", "b(sympy_lower,sympy_upper,1)"),
+     ("b(::)", "b(sympy_lower,sympy_upper,1)"),
+     ("b(5::)", "b(5,sympy_upper,1)"),
+     ("b(:5:)", "b(sympy_lower,5,1)"),
+     ("b(::5)", "b(sympy_lower,sympy_upper,5)"),
+     ("b(i::)", "b(i,sympy_upper,1)"),
+     ("b(:i:)", "b(sympy_lower,i,1)"),
+     ("b(::i)", "b(sympy_lower,sympy_upper,i)"),
+     ("b(i:5:)", "b(i,5,1)"),
+     ("b(i:j:)", "b(i,j,1)"),
+     ("b(i::j)", "b(i,sympy_upper,j)"),
+     ("b(:i:j)", "b(sympy_lower,i,j)"),
+     ("b(i:j:k)", "b(i,j,k)"),
+     ("b", "b(sympy_no_bounds,sympy_no_bounds,1)"),
+     ("c(i,j)", "c(i,i,1,j,j,1)"),
+     ("c(::,::)", "c(sympy_lower,sympy_upper,1,"
+      "sympy_lower,sympy_upper,1)"),
+     ("c", "c(sympy_no_bounds,sympy_no_bounds,1,sympy_no_bounds,"
+      "sympy_no_bounds,1)"),
+     ("b(i)%x", "b_x(i,i,1,sympy_no_bounds,sympy_no_bounds,1)"),
+     ("b(i)%x(j)", "b_x(i,i,1,j,j,1)"),
+     ("c(i,j)%x", "c_x(i,i,1,j,j,1,sympy_no_bounds,sympy_no_bounds,1)"),
+     ("c(i,j)%x(j)", "c_x(i,i,1,j,j,1,j,j,1)"),
+     ("c(i,j)%d%e", "c_d_e(i,i,1,j,j,1,sympy_no_bounds,sympy_no_bounds,1,"
+      "sympy_no_bounds,sympy_no_bounds,1)"),
+     ("c(i,j)%d%f(i)", "c_d_f(i,i,1,j,j,1,sympy_no_bounds,sympy_no_bounds,1,"
+      "i,i,1)"),
+     ("c(i::k,j)%d%f(i:j:k)",
+      "c_d_f(i,sympy_upper,k,j,j,1,sympy_no_bounds,sympy_no_bounds,1,i,j,k)"),
+     # Check name clashes, if a user variable is the same as the names
+     # for upper/lower bound
+     ("sympy_upper(:)", "sympy_upper(sympy_lower,sympy_upper_1,1)"),
+     ("sympy_lower(:)", "sympy_lower(sympy_lower_1,sympy_upper,1)"),
+     # The +sympy_upper at the end is an array expression but should not get
+     # indices added!
+     ("sympy_lower(:)+sympy_upper",
+      "sympy_lower(sympy_lower_1,sympy_upper_1,1) + sympy_upper"
+      "(sympy_no_bounds,sympy_no_bounds,1)"),
+     ])
 def test_sym_writer_array_expressions(fortran_reader, expressions):
     '''Test that array expressions (including ones using user-defined
     types) are converted correctly. A Fortran range is converted into
@@ -465,20 +459,24 @@ def test_gen_indices():
     assert "unsupported gen_indices index 'None'" in str(err.value)
 
 
-@pytest.mark.parametrize("fortran_expr,sympy_str",
-                         [("a%b", "a_b"),
-                          # Handle name clash:
-                          ("a%c + a_c", "a_c_1 + a_c"),
-                          ("a%b(i)", "a_b(i,i,1)"),
-                          ("b(i)%b", "b_b(i,i,1)"),
-                          ("b(:)%b(i) + b(1)%c",
-                           "b_b(sympy_lower,sympy_upper,1,i,i,1) + "
-                           "b_c(1,1,1)"),
-                          ("b(i)%b(j)", "b_b(i,i,1,j,j,1)"),
-                          ("a%b(i)%c", "a_b_c(i,i,1)"),
-                          ("a%b%c(i)", "a_b_c(i,i,1)"),
-                          ("a%b%c(2 * i - 1)", "a_b_c(2 * i - 1,2 * i - 1,1)")
-                          ])
+@pytest.mark.parametrize(
+    "fortran_expr,sympy_str",
+    [("a%b", "a_b"),
+     # Handle name clash:
+     ("a%c + a_c", "a_c + a_c_1"),
+     ("a%b(i)", "a_b(sympy_no_bounds,sympy_no_bounds,1,i,i,1)"),
+     ("b(i)%b", "b_b(i,i,1,sympy_no_bounds,sympy_no_bounds,1)"),
+     ("b(:)%b(i) + b(1)%c",
+      "b_b(sympy_lower,sympy_upper,1,i,i,1) + "
+      "b_c(1,1,1,sympy_no_bounds,sympy_no_bounds,1)"),
+     ("b(i)%b(j)", "b_b(i,i,1,j,j,1)"),
+     ("a%b(i)%c", "a_b_c(sympy_no_bounds,sympy_no_bounds,1,i,i,1,"
+      "sympy_no_bounds,sympy_no_bounds,1)"),
+     ("a%b%c(i)", "a_b_c(sympy_no_bounds,sympy_no_bounds,1,sympy_no_bounds,"
+      "sympy_no_bounds,1,i,i,1)"),
+     ("a%b%c(2 * i - 1)", "a_b_c(sympy_no_bounds,sympy_no_bounds,1,"
+      "sympy_no_bounds,sympy_no_bounds,1,2 * i - 1,2 * i - 1,1)")
+     ])
 def test_sympy_writer_user_types(fortran_reader, fortran_writer,
                                  fortran_expr, sympy_str):
     '''Test handling of user-defined types, e.g. conversion of
diff --git a/src/psyclone/tests/psyir/frontend/fparser2_where_handler_test.py b/src/psyclone/tests/psyir/frontend/fparser2_where_handler_test.py
index 866dccb69..894ee8120 100644
--- a/src/psyclone/tests/psyir/frontend/fparser2_where_handler_test.py
+++ b/src/psyclone/tests/psyir/frontend/fparser2_where_handler_test.py
@@ -1177,8 +1177,8 @@ def test_elemental_function_to_loop(fortran_reader, fortran_writer):
     psyir = fortran_reader.psyir_from_source(code)
     assert isinstance(psyir.children[0].children[0], CodeBlock)
     correct = '''! PSyclone CodeBlock (unsupported code) reason:
-  !  - PSyclone doesn't yet support reference to imported \
-symbols inside WHERE clauses.
+  !  - PSyclone doesn't yet support references to imported/unresolved \
+symbols inside WHERE clauses: 'somefunc' is unresolved.
   WHERE (somefunc(a) < 2)
     b = a
   END WHERE'''
@@ -1196,8 +1196,8 @@ symbols inside WHERE clauses.
     psyir = fortran_reader.psyir_from_source(code)
     assert isinstance(psyir.children[0].children[0], CodeBlock)
     correct = '''! PSyclone CodeBlock (unsupported code) reason:
-  !  - PSyclone doesn't yet support reference to imported \
-symbols inside WHERE clauses.
+  !  - PSyclone doesn't yet support references to imported/unresolved \
+symbols inside WHERE clauses: 'somefunc' is unresolved.
   WHERE (somefunc(a) < 2)
     b = a
   END WHERE'''
diff --git a/src/psyclone/tests/psyir/frontend/sympy_reader_test.py b/src/psyclone/tests/psyir/frontend/sympy_reader_test.py
index 539686b51..a792d89f2 100644
--- a/src/psyclone/tests/psyir/frontend/sympy_reader_test.py
+++ b/src/psyclone/tests/psyir/frontend/sympy_reader_test.py
@@ -63,7 +63,7 @@ def test_sympy_reader_constructor():
                                          ("c(i,j)", "c(i,j)"),
                                          ("b(2:3:4)", "b(2:3:4)"),
                                          ("b(2:3:1)", "b(2:3)"),
-                                         ("b", "b(:)"),
+                                         ("b", "b"),
                                          ("b(:)", "b(:)"),
                                          ("b(::)", "b(:)"),
                                          ("b(::1)", "b(:)"),
diff --git a/src/psyclone/tests/psyir/transformations/intrinsics/minval2loop_trans_test.py b/src/psyclone/tests/psyir/transformations/intrinsics/minval2loop_trans_test.py
index d644b6817..bf0225d99 100644
--- a/src/psyclone/tests/psyir/transformations/intrinsics/minval2loop_trans_test.py
+++ b/src/psyclone/tests/psyir/transformations/intrinsics/minval2loop_trans_test.py
@@ -68,9 +68,7 @@ def test_init_var():
     trans = Minval2LoopTrans()
     var_symbol = DataSymbol("var", REAL_TYPE)
     result = trans._init_var(Reference(var_symbol))
-    # As 'huge' is not yet part of an expression, the 'debug_string()'
-    # method incorrectly assumes it is a call.
-    assert result.debug_string() == "call HUGE(var)\n"
+    assert result.debug_string() == "HUGE(var)\n"
 
 
 def test_str():
diff --git a/src/psyclone/tests/psyir/transformations/reference2arrayrange_trans_test.py b/src/psyclone/tests/psyir/transformations/reference2arrayrange_trans_test.py
index 90b12b423..a006a5bde 100644
--- a/src/psyclone/tests/psyir/transformations/reference2arrayrange_trans_test.py
+++ b/src/psyclone/tests/psyir/transformations/reference2arrayrange_trans_test.py
@@ -194,18 +194,19 @@ def test_multid(fortran_reader, fortran_writer):
 
 
 def test_intrinsics(fortran_reader, fortran_writer):
-    '''Test that references to arrays within intrinsics are transformed to
+    '''Test that references to arrays within intrinsics are not transformed to
     array slice notation, using dotproduct as the example.
 
     '''
     code = CODE.replace("a = b", "b = dot_product(a, a(:))")
     result = apply_trans(fortran_reader, fortran_writer, code)
-    assert "b = DOT_PRODUCT(a(:), a(:))" in result
+    assert "b = DOT_PRODUCT(a, a(:))" in result
 
 
 def test_call(fortran_reader, fortran_writer):
     '''Test that references to arrays that are arguments to a call are
-    transformed to array slice notation.
+    *not* transformed to array slice notation (since this affects the
+    bounds of the array seen within the called routine).
 
     '''
     code = (
@@ -216,7 +217,28 @@ def test_call(fortran_reader, fortran_writer):
         "  call work(a,b)\n"
         "end program test\n")
     result = apply_trans(fortran_reader, fortran_writer, code)
-    assert "call work(a(:), b)" in result
+    assert "call work(a, b)" in result
+
+
+def test_ambiguous_call_array_reference(fortran_reader, fortran_writer):
+    '''
+    Test that references to arrays that *may* be arguments to a call are *not*
+    transformed because, if it is a call, this will affect the bounds of the
+    array seen within the called routine.
+
+    '''
+    code = '''\
+    program test
+      use some_mod, only: work
+      integer, dimension(10) :: a, b
+      a = 1
+      ! Without resolving 'work', we don't know whether it is an array
+      ! or a function.
+      b = work(a)
+    end program test'''
+    result = apply_trans(fortran_reader, fortran_writer, code)
+    assert "a(:) = 1" in result
+    assert "b(:) = work(a)" in result
 
 
 def test_validate():
@@ -273,9 +295,9 @@ def test_validate_query(fortran_reader):
         for reference in location.walk(Reference)[1:]:
             with pytest.raises(TransformationError) as info:
                 trans.validate(reference)
-            assert (f"References to arrays passed as arguments to intrinsic "
-                    f"enquiry routine '{text}' should not be transformed."
-                    in str(info.value))
+            assert (f"supplied node is passed as an argument to a Call to a "
+                    f"non-elemental routine ({text}(a, 1)) and should not be "
+                    f"transformed." in str(info.value))
 
     # Check the references to 'b' in the hidden lbound and ubound
     # intrinsics within 'b(:)' do not get modified.
@@ -287,24 +309,25 @@ def test_validate_query(fortran_reader):
                 isinstance(reference.symbol, IntrinsicSymbol)):
             with pytest.raises(TransformationError) as info:
                 trans.validate(reference)
-            assert ("References to arrays passed as arguments to intrinsic "
-                    "enquiry routine '" in str(info.value))
+            assert ("supplied node is passed as an argument to a Call to a "
+                    "non-elemental routine (" in str(info.value))
 
     # Check the reference to 'b' in the size intrinsics does not get modified
     assignment = psyir.children[0].children[2]
     reference = assignment.children[1].arguments[0]
     with pytest.raises(TransformationError) as info:
         trans.validate(reference)
-    assert ("References to arrays passed as arguments to intrinsic enquiry "
-            "routine 'SIZE' should not be transformed." in str(info.value))
+    assert ("supplied node is passed as an argument to a Call to a "
+            "non-elemental routine (SIZE" in str(info.value))
 
     ifblock = psyir.walk(IfBlock)[0]
     allocd = ifblock.condition
     assert isinstance(allocd, IntrinsicCall)
     with pytest.raises(TransformationError) as info:
         trans.validate(allocd.arguments[0])
-    assert ("References to arrays passed as arguments to intrinsic enquiry "
-            "routine 'ALLOCATED' should not be transformed" in str(info.value))
+    assert ("supplied node is passed as an argument to a Call to a "
+            "non-elemental routine (ALLOCATED(igor)) and should not be "
+            "transformed" in str(info.value))
 
 
 def test_validate_structure(fortran_reader):
@@ -345,8 +368,8 @@ def test_validate_deallocate(fortran_reader):
     trans = Reference2ArrayRangeTrans()
     with pytest.raises(TransformationError) as info:
         trans.validate(reference)
-    assert ("References to arrays passed to 'DEALLOCATE' intrinsics should not"
-            " be transformed, but found:\n DEALLOCATE(a)" in str(info.value))
+    assert ("node is passed as an argument to a Call to a non-elemental "
+            "routine (DEALLOCATE(a)" in str(info.value))
 
 
 def test_validate_pointer_assignment(fortran_reader):