Repository: stfc/PSyclone
Pull request #2372 (2023-10-19 16:00:56)

(Towards #1960) put WHEREs containing reductions into CodeBlocks and fix non-unity lower-bound bug.


Patch:
diff --git a/changelog b/changelog
index d59f6ba73..95be5eaca 100644
--- a/changelog
+++ b/changelog
@@ -152,6 +152,8 @@
 	55) PR #2465 towards #2004. Move same_range functionality to
 	ArrayMixin and other improvements.
 
+	56) PR #2372 towards #1960. Fix bugs on Fortran WHERE canonicalisation. 
+
 release 2.4.0 29th of September 2023
 
 	1) PR #1758 for #1741. Splits the PSyData read functionality into a
diff --git a/examples/nemo/scripts/process_nemo.py b/examples/nemo/scripts/process_nemo.py
index 841bb47e6..c73c996ad 100755
--- a/examples/nemo/scripts/process_nemo.py
+++ b/examples/nemo/scripts/process_nemo.py
@@ -78,10 +78,8 @@ NOT_PERFORMANT = [
 # Files that we won't touch at all, either because PSyclone actually fails
 # or because it produces incorrect Fortran.
 NOT_WORKING = [
-    # Array accessed inside WHERE does not use array notation
+    # TODO #717 - array accessed inside WHERE does not use array notation
     "diurnal_bulk.f90",
-    # TODO #1960: There is a bug on a WHERE fparser2 frontend processing
-    "sbccpl.f90",
 ]
 
 
diff --git a/src/psyclone/psyir/frontend/fparser2.py b/src/psyclone/psyir/frontend/fparser2.py
index e6afda7d5..2e30c32c3 100644
--- a/src/psyclone/psyir/frontend/fparser2.py
+++ b/src/psyclone/psyir/frontend/fparser2.py
@@ -3604,15 +3604,17 @@ class Fparser2Reader():
         '''
         Utility function that modifies each ArrayReference object in the
         supplied PSyIR fragment so that they are indexed using the supplied
-        loop variables rather than having colon array notation.
+        loop variables rather than having colon array notation. This indexing
+        is always done relative to the declared lower bound of the array being
+        accessed.
 
-        :param parent: root of PSyIR sub-tree to search for Array \
+        :param parent: root of PSyIR sub-tree to search for Array
                        references to modify.
         :type parent:  :py:class:`psyclone.psyir.nodes.Node`
         :param loop_vars: the variable names for the array indices.
         :type loop_vars: list of str
 
-        :raises NotImplementedError: if array sections of differing ranks are \
+        :raises NotImplementedError: if array sections of differing ranks are
                                      found.
         '''
         assigns = parent.walk(Assignment)
@@ -3629,6 +3631,8 @@ class Fparser2Reader():
         # PSyIR (using e.g. the Fortran backend) will not
         # compile. We need to implement robust identification of the
         # types of all symbols in the PSyIR fragment.
+        table = parent.scope.symbol_table
+        one = Literal("1", INTEGER_TYPE)
         arrays = parent.walk(ArrayMixin)
         first_rank = None
         for array in arrays:
@@ -3650,15 +3654,58 @@ class Fparser2Reader():
             else:
                 first_rank = rank
 
-            # Replace the PSyIR Ranges with the loop variables
+            base_ref = _copy_full_base_reference(array)
+            array_ref = array.ancestor(Reference, include_self=True)
+            shape = array_ref.datatype.shape
+            add_op = BinaryOperation.Operator.ADD
+            sub_op = BinaryOperation.Operator.SUB
+            # Replace the PSyIR Ranges with appropriate index expressions.
             range_idx = 0
             for idx, child in enumerate(array.indices):
-                if isinstance(child, Range):
-                    symbol = _find_or_create_unresolved_symbol(
-                        array, loop_vars[range_idx],
-                        symbol_type=DataSymbol, datatype=UnresolvedType())
-                    array.children[idx] = Reference(symbol)
-                    range_idx += 1
+                if not isinstance(child, Range):
+                    continue
+                # We need the lower bound of the appropriate dimension of this
+                # array as we will index relative to it. Note that the 'shape'
+                # of the datatype only gives us extents, not the lower bounds
+                # of the declaration or slice.
+                if isinstance(shape[range_idx], ArrayType.Extent):
+                    # We don't know the bounds of this array so we have
+                    # to query using LBOUND.
+                    lbound = IntrinsicCall.create(
+                        IntrinsicCall.Intrinsic.LBOUND,
+                        [base_ref.copy(),
+                         ("dim", Literal(str(idx+1), INTEGER_TYPE))])
+                else:
+                    if array.is_full_range(idx):
+                        # The access to this index is to the full range of
+                        # the array.
+                        # TODO #949 - ideally we would try to find the lower
+                        # bound of the array by interrogating `array.symbol.
+                        # datatype` but the fparser2 frontend doesn't currently
+                        # support array declarations with explicit lower bounds
+                        lbound = IntrinsicCall.create(
+                            IntrinsicCall.Intrinsic.LBOUND,
+                            [base_ref.copy(),
+                             ("dim", Literal(str(idx+1), INTEGER_TYPE))])
+                    else:
+                        # We need the lower bound of this access.
+                        lbound = child.start.copy()
+
+                # Create the index expression.
+                symbol = table.lookup(loop_vars[range_idx])
+                if isinstance(lbound, Literal) and lbound.value == "1":
+                    # Lower bound is just unity so we can use the loop-idx
+                    # directly.
+                    expr2 = Reference(symbol)
+                else:
+                    # We don't know what the lower bound is so have to
+                    # have an expression:
+                    #    idx-expr = array-lower-bound + loop-idx - 1
+                    expr = BinaryOperation.create(
+                        add_op, lbound, Reference(symbol))
+                    expr2 = BinaryOperation.create(sub_op, expr, one.copy())
+                array.children[idx] = expr2
+                range_idx += 1
 
     def _where_construct_handler(self, node, parent):
         '''
@@ -3678,7 +3725,7 @@ class Fparser2Reader():
             WHERE(logical-mask) statement
 
         :param node: node in the fparser2 parse tree representing the WHERE.
-        :type node: :py:class:`fparser.two.Fortran2003.Where_Construct` or \
+        :type node: :py:class:`fparser.two.Fortran2003.Where_Construct` |
                     :py:class:`fparser.two.Fortran2003.Where_Stmt`
         :param parent: parent node in the PSyIR.
         :type parent: :py:class:`psyclone.psyir.nodes.Node`
@@ -3686,11 +3733,42 @@ class Fparser2Reader():
         :returns: the top-level Loop object in the created loop nest.
         :rtype: :py:class:`psyclone.psyir.nodes.Loop`
 
-        :raises InternalError: if the parse tree does not have the expected \
-                               structure.
-        :raises NotImplementedError: if the logical mask of the WHERE does \
-                                     not use array notation.
+        :raises InternalError: if the parse tree does not have the expected
+            structure.
+        :raises NotImplementedError: if the parse tree contains a Fortran
+            intrinsic that performs a reduction but still returns an array.
+        :raises NotImplementedError: if the logical mask of the WHERE does
+            not use array notation.
+
         '''
+        def _contains_intrinsic_reduction(pnodes):
+            '''
+            Utility to check for Fortran intrinsics that perform a reduction
+            but result in an array.
+
+            :param pnodes: node(s) in the parse tree to check.
+            :type pnodes: list[:py:class:`fparser.two.utils.Base`] |
+                          :py:class:`fparser.two.utils.Base`
+
+            :returns: whether or not the supplied node(s) in the parse tree
+                      contain a call to an intrinsic that performs a reduction
+                      into an array.
+            :rtype: bool
+
+            '''
+            intr_nodes = walk(pnodes, Fortran2003.Intrinsic_Function_Reference)
+            for intr in intr_nodes:
+                if (intr.children[0].string in
+                        Fortran2003.Intrinsic_Name.array_reduction_names):
+                    # These intrinsics are only a problem if they return an
+                    # array rather than a scalar.
+                    arg_specs = walk(intr.children[1],
+                                     Fortran2003.Actual_Arg_Spec)
+                    if any(spec.children[0].string == 'dim'
+                           for spec in arg_specs):
+                        return True
+            return False
+
         if isinstance(node, Fortran2003.Where_Stmt):
             # We have a Where statement. Check that the parse tree has the
             # expected structure.
@@ -3704,6 +3782,11 @@ class Fparser2Reader():
                     f"Expected the second entry of a Fortran2003.Where_Stmt "
                     f"items tuple to be an Assignment_Stmt but found: "
                     f"{type(node.items[1]).__name__}")
+            if _contains_intrinsic_reduction(node.items[1]):
+                raise NotImplementedError(
+                    f"TODO #1960 - WHERE statements which contain array-"
+                    f"reduction intrinsics are not supported but found "
+                    f"'{node}'")
             was_single_stmt = True
             annotations = ["was_where", "was_single_stmt"]
             logical_expr = [node.items[0]]
@@ -3717,6 +3800,11 @@ class Fparser2Reader():
             if not isinstance(node.content[-1], Fortran2003.End_Where_Stmt):
                 raise InternalError(f"Failed to find closing end where "
                                     f"statement in: {node}")
+            if _contains_intrinsic_reduction(node.content[1:-1]):
+                raise NotImplementedError(
+                    f"TODO #1960 - WHERE constructs which contain an array-"
+                    f"reduction intrinsic are not supported but found "
+                    f"'{node}'")
             was_single_stmt = False
             annotations = ["was_where"]
             logical_expr = node.content[0].items
@@ -3729,6 +3817,12 @@ class Fparser2Reader():
         # the NEMO style where the fact that `a` is an array is made
         # explicit using the colon notation, e.g. `a(:, :) < 0.0`.
 
+        if _contains_intrinsic_reduction(logical_expr):
+            raise NotImplementedError(
+                f"TODO #1960 - WHERE constructs which contain an array-"
+                f"reduction intrinsic in their logical expression are not "
+                f"supported but found '{logical_expr}'")
+
         # For this initial processing of the logical-array expression we
         # use a temporary parent as we haven't yet constructed the PSyIR
         # for the loop nest and innermost IfBlock. Once we have a valid
@@ -3742,10 +3836,11 @@ class Fparser2Reader():
             # because the code doesn't use explicit array syntax. At least one
             # variable in the logical-array expression must be an array for
             # this to be a valid WHERE().
-            # TODO #717. Look-up the shape of the array in the SymbolTable.
+            # TODO #1799. Look-up the shape of the array in the SymbolTable.
             raise NotImplementedError(
                 f"Only WHERE constructs using explicit array notation (e.g. "
                 f"my_array(:,:)) are supported but found '{logical_expr}'.")
+
         for array in arrays:
             if any(isinstance(idx, Range) for idx in array.indices):
                 first_array = array
@@ -3753,12 +3848,14 @@ class Fparser2Reader():
         else:
             raise NotImplementedError(
                 f"Only WHERE constructs using explicit array notation "
-                f"including ranges (e.g. my_array(1,:) are supported but "
+                f"including ranges (e.g. 'my_array(1,:)') are supported but "
                 f"found '{logical_expr}'")
 
-        # All array sections in a Fortran WHERE must have the same rank so
-        # just look at the first array.
-        rank = self._array_notation_rank(first_array)
+        array_ref = first_array.ancestor(Reference, include_self=True)
+        mask_shape = array_ref.datatype.shape
+        # All array sections in a Fortran WHERE must have the same shape so
+        # just look at that of the mask.
+        rank = len(mask_shape)
         # Create a list to hold the names of the loop variables as we'll
         # need them to index into the arrays.
         loop_vars = rank*[""]
@@ -3778,35 +3875,23 @@ class Fparser2Reader():
                         annotations=annotations)
             # Point to the original WHERE statement in the parse tree.
             loop.ast = node
-            # Add loop lower bound
+
+            # This loop is over the *shape* of the mask and thus starts
+            # at unity. Each individual array access is then adjusted
+            # according to the lower bound of that array.
             loop.addchild(Literal("1", integer_type))
-            # Add loop upper bound - we use the SIZE operator to query the
-            # extent of the current array dimension
-            size_node = IntrinsicCall(IntrinsicCall.Intrinsic.SIZE,
-                                      parent=loop)
-            loop.addchild(size_node)
-
-            # Create the first argument to the SIZE operator
-            if isinstance(first_array, Member):
-                # The array access is a member of some derived type
-                parent_ref = first_array.ancestor(Reference)
-                new_ref = parent_ref.copy()
-                orig_member = parent_ref.member
-                member = new_ref.member
-                while orig_member is not first_array:
-                    member = member.member
-                    orig_member = orig_member.member
-                member.parent.children[0] = Member(first_array.name,
-                                                   parent=member.parent)
+            # Add loop upper bound using the shape of the mask.
+            if isinstance(mask_shape[idx-1], ArrayType.Extent):
+                # We don't have an explicit value for the upper bound so we
+                # have to query it using SIZE.
+                loop.addchild(
+                    IntrinsicCall.create(IntrinsicCall.Intrinsic.SIZE,
+                                         [array_ref.copy(),
+                                          ("dim", Literal(str(idx),
+                                                          integer_type))]))
             else:
-                # The array access is to a symbol of ArrayType
-                symbol = _find_or_create_unresolved_symbol(
-                    size_node, first_array.name, symbol_type=DataSymbol,
-                    datatype=UnresolvedType())
-                new_ref = Reference(symbol)
-            size_node.addchild(new_ref)
-            size_node.addchild(Literal(str(idx), integer_type,
-                                       parent=size_node))
+                loop.addchild(mask_shape[idx-1].upper.copy())
+
             # Add loop increment
             loop.addchild(Literal("1", integer_type))
             # Fourth child of a Loop must be a Schedule
@@ -3822,6 +3907,7 @@ class Fparser2Reader():
                 # handler returns
                 root_loop = loop
             new_parent = sched
+
         # Now we have the loop nest, add an IF block to the innermost
         # schedule
         ifblock = IfBlock(parent=new_parent, annotations=annotations)
diff --git a/src/psyclone/psyir/nodes/array_mixin.py b/src/psyclone/psyir/nodes/array_mixin.py
index d1ccab239..640f03bcc 100644
--- a/src/psyclone/psyir/nodes/array_mixin.py
+++ b/src/psyclone/psyir/nodes/array_mixin.py
@@ -516,48 +516,70 @@ class ArrayMixin(metaclass=abc.ABCMeta):
                     f"expression but found '{type(child).__name__}'")
         return self.children
 
+    def _extent(self, idx):
+        '''
+        Create PSyIR for the number of elements in dimension `idx` of this
+        array access. It is given by (stop - start)/step + 1 or, if it is for
+        the full range, by the SIZE intrinsic.
+
+        :param int idx: the array index for which to compute the number of
+                        elements.
+
+        :returns: the PSyIR expression for the number of elements in the
+                  specified array index.
+        :rtype: :py:class:`psyclone.psyir.nodes.BinaryOperation` |
+                :py:class:`psyclone.psyir.nodes.IntrinsicCall`
+        '''
+        expr = self.indices[idx]
+        one = Literal("1", INTEGER_TYPE)
+
+        if isinstance(expr, Range):
+            start = expr.start
+            stop = expr.stop
+            step = expr.step
+        else:
+            # No range so just a single element is accessed.
+            return one
+
+        if (isinstance(start, IntrinsicCall) and
+                isinstance(stop, IntrinsicCall) and self.is_full_range(idx)):
+            # Access is to full range and start and stop are expressed in terms
+            # of LBOUND and UBOUND. Therefore, it's simpler to use SIZE.
+            return IntrinsicCall.create(
+                IntrinsicCall.Intrinsic.SIZE,
+                [start.children[0].copy(),
+                 ("dim", Literal(str(idx+1), INTEGER_TYPE))])
+
+        if start == one and step == one:
+            # The range starts at 1 and the step is 1 so the extent is just
+            # the upper bound.
+            return stop.copy()
+
+        extent = BinaryOperation.create(BinaryOperation.Operator.SUB,
+                                        stop.copy(), start.copy())
+        if step != one:
+            # Step is not unity so have to divide range by it.
+            result = BinaryOperation.create(BinaryOperation.Operator.DIV,
+                                            extent, step.copy())
+        else:
+            result = extent
+        # Extent is currently 'stop-start' or '(stop-start)/step' so we have
+        # to add a '+ 1'
+        return BinaryOperation.create(BinaryOperation.Operator.ADD,
+                                      result, one.copy())
+
     def _get_effective_shape(self):
         '''
         :returns: the shape of the array access represented by this node.
-        :rtype: List[:py:class:`psyclone.psyir.nodes.DataNode`]
+        :rtype: list[:py:class:`psyclone.psyir.nodes.DataNode`]
 
         :raises NotImplementedError: if any of the array-indices involve a
                                      function call or an expression.
         '''
-        def _num_elements(expr):
-            '''
-            Create PSyIR for the number of elements in this range. It
-            is given by (stop - start)/step + 1.
-
-            :param expr: the range for which to compute the number of elements.
-            :type expr: :py:class:`psyclone.psyir.nodes.Range` or \
-                :py:class:`psyclone.psyir.symbols.ArrayType.ArrayBounds`
-
-            :returns: the PSyIR expression for the number of elements in the \
-                      supplied range.
-            :rtype: :py:class:`psyclone.psyir.nodes.BinaryOperation`
-
-            '''
-            if isinstance(expr, Range):
-                start = expr.start
-                stop = expr.stop
-                step = expr.step
-            elif isinstance(expr, ArrayType.ArrayBounds):
-                start = expr.lower
-                stop = expr.upper
-                step = Literal("1", INTEGER_TYPE)
-            minus = BinaryOperation.create(BinaryOperation.Operator.SUB,
-                                           stop.copy(), start.copy())
-            div = BinaryOperation.create(BinaryOperation.Operator.DIV,
-                                         minus, step.copy())
-            plus = BinaryOperation.create(BinaryOperation.Operator.ADD,
-                                          div, Literal("1", INTEGER_TYPE))
-            return plus
-
         shape = []
-        for idx_expr in self.indices:
+        for idx, idx_expr in enumerate(self.indices):
             if isinstance(idx_expr, Range):
-                shape.append(_num_elements(idx_expr))
+                shape.append(self._extent(idx))
 
             elif isinstance(idx_expr, Reference):
                 dtype = idx_expr.datatype
@@ -571,7 +593,9 @@ class ArrayMixin(metaclass=abc.ABCMeta):
                             f"another array must be 1D but '{idx_expr.name}' "
                             f"used to index into '{self.name}' has "
                             f"{len(indirect_array_shape)} dimensions.")
-                    shape.append(_num_elements(dtype.shape[0]))
+                    # pylint: disable=protected-access
+                    shape.append(idx_expr._extent(idx))
+
             elif isinstance(idx_expr, (Call, Operation, CodeBlock)):
                 # We can't yet straightforwardly query the type of a function
                 # call or Operation - TODO #1799.
diff --git a/src/psyclone/psyir/nodes/array_reference.py b/src/psyclone/psyir/nodes/array_reference.py
index 4d815e5e1..11c6461b0 100644
--- a/src/psyclone/psyir/nodes/array_reference.py
+++ b/src/psyclone/psyir/nodes/array_reference.py
@@ -47,7 +47,7 @@ from psyclone.psyir.nodes.reference import Reference
 from psyclone.psyir.symbols import (DataSymbol, UnresolvedType,
                                     UnsupportedFortranType, UnsupportedType,
                                     DataTypeSymbol, ScalarType, ArrayType,
-                                    INTEGER_TYPE)
+                                    INTEGER_TYPE, Symbol)
 
 
 class ArrayReference(ArrayMixin, Reference):
@@ -87,7 +87,7 @@ class ArrayReference(ArrayMixin, Reference):
                 f"indices argument in create method of ArrayReference class "
                 f"should be a list but found '{type(indices).__name__}'.")
         if not symbol.is_array:
-            # Deferred and Unsupported types may still be arrays
+            # Unresolved and Unsupported types may still be arrays
             if not isinstance(symbol.datatype, (UnresolvedType,
                                                 UnsupportedType)):
                 raise GenerationError(
@@ -130,7 +130,11 @@ class ArrayReference(ArrayMixin, Reference):
         '''
         shape = self._get_effective_shape()
         if shape:
-            if isinstance(self.symbol.datatype, ArrayType):
+            if type(self.symbol) is Symbol:
+                # We don't have any information on the shape of the original
+                # declaration.
+                orig_shape = None
+            elif isinstance(self.symbol.datatype, ArrayType):
                 # We have full type information so we know the shape of the
                 # original declaration.
                 orig_shape = self.symbol.datatype.shape
@@ -141,9 +145,9 @@ class ArrayReference(ArrayMixin, Reference):
                 orig_shape = self.symbol.datatype.partial_datatype.shape
             else:
                 # We don't have any information on the shape of the original
-                # delcaration.
-                orig_shape = []
-            if (len(shape) == len(orig_shape) and
+                # declaration.
+                orig_shape = None
+            if (orig_shape is not None and len(shape) == len(orig_shape) and
                     all(self.is_full_range(idx) for idx in range(len(shape)))):
                 # Although this access has a shape, it is in fact for the
                 # whole array and therefore the type of the result is just
@@ -151,8 +155,9 @@ class ArrayReference(ArrayMixin, Reference):
                 # StructureReference but they have their own implementation
                 # of this method.)
                 return self.symbol.datatype
-            if isinstance(self.symbol.datatype, UnsupportedType):
-                # Even if an Unknown(Fortran)Type has partial type
+            if type(self.symbol) is Symbol or isinstance(self.symbol.datatype,
+                                                         UnsupportedType):
+                # Even if an Unsupported(Fortran)Type has partial type
                 # information, we can't easily use it here because we'd need
                 # to re-write the original Fortran declaration stored in the
                 # type. We could manipulate the shape in the fparser2 parse
@@ -167,6 +172,8 @@ class ArrayReference(ArrayMixin, Reference):
             return ArrayType(base_type, shape)
 
         # Otherwise, we're accessing a single element of the array.
+        if type(self.symbol) is Symbol:
+            return UnresolvedType()
         if isinstance(self.symbol.datatype, UnsupportedType):
             if (isinstance(self.symbol.datatype, UnsupportedFortranType) and
                     self.symbol.datatype.partial_datatype):
diff --git a/src/psyclone/psyir/nodes/structure_reference.py b/src/psyclone/psyir/nodes/structure_reference.py
index 30e0d4d70..f3c920c81 100644
--- a/src/psyclone/psyir/nodes/structure_reference.py
+++ b/src/psyclone/psyir/nodes/structure_reference.py
@@ -303,7 +303,7 @@ class StructureReference(Reference):
         if isinstance(dtype, (UnresolvedType, UnsupportedType)):
             # We don't know the type of the symbol that defines the type
             # of this structure.
-            return UnresolvedType()
+            dtype = UnresolvedType()
 
         # We do have the definition of this structure - walk down it.
         cursor = self
@@ -325,9 +325,10 @@ class StructureReference(Reference):
                 cursor_type = cursor_type.intrinsic
             if isinstance(cursor_type, DataTypeSymbol):
                 cursor_type = cursor_type.datatype
-            cursor_type = cursor_type.components[cursor.name].datatype
-            if isinstance(cursor_type, (UnsupportedType, UnresolvedType)):
-                return UnresolvedType()
+            if not isinstance(cursor_type, (UnresolvedType, UnsupportedType)):
+                # Once we've hit an Unresolved/UnsupportedType the cursor_type
+                # will remain set to that as we can't do any better.
+                cursor_type = cursor_type.components[cursor.name].datatype
             if isinstance(cursor, ArrayMixin):
                 # pylint: disable=protected-access
                 shape.extend(cursor._get_effective_shape())
@@ -346,7 +347,7 @@ class StructureReference(Reference):
                 else:
                     # No indices so it is an access to a whole array.
                     cursor_shape = cursor_type.shape
-                if cursor_shape and shape != cursor_shape:
+                if cursor_shape and len(shape) != len(cursor_shape):
                     # This ultimate access is an array but we've already
                     # encountered one or more slices earlier in the access
                     # expression.
diff --git a/src/psyclone/psyir/transformations/inline_trans.py b/src/psyclone/psyir/transformations/inline_trans.py
index 12abe09d9..bf27c9b4d 100644
--- a/src/psyclone/psyir/transformations/inline_trans.py
+++ b/src/psyclone/psyir/transformations/inline_trans.py
@@ -776,8 +776,11 @@ class InlineTrans(Transformation):
             # TODO #924. It would be useful if the `datatype` property was
             # a method that took an optional 'resolve' argument to indicate
             # that it should attempt to resolve any UnresolvedTypes.
-            if isinstance(actual_arg.datatype, (UnresolvedType,
-                                                UnsupportedType)):
+            if (isinstance(actual_arg.datatype,
+                           (UnresolvedType, UnsupportedType)) or
+                (isinstance(actual_arg.datatype, ArrayType) and
+                 isinstance(actual_arg.datatype.intrinsic,
+                            (UnresolvedType, UnsupportedType)))):
                 raise TransformationError(
                     f"Routine '{routine.name}' cannot be inlined because "
                     f"the type of the actual argument "


Test patch:
diff --git a/src/psyclone/tests/nemo/transformations/openacc/kernels_directive_test.py b/src/psyclone/tests/nemo/transformations/openacc/kernels_directive_test.py
index b8f00d930..ed7f84bfa 100644
--- a/src/psyclone/tests/nemo/transformations/openacc/kernels_directive_test.py
+++ b/src/psyclone/tests/nemo/transformations/openacc/kernels_directive_test.py
@@ -254,10 +254,12 @@ def test_kernels_around_where_construct(parser):
     assert isinstance(schedule[0].dir_body[0], Loop)
     new_code = str(psy.gen)
     assert ("  !$acc kernels\n"
-            "  do widx2 = 1, SIZE(a, 2), 1\n"
-            "    do widx1 = 1, SIZE(a, 1), 1\n"
-            "      if (a(widx1,widx2) < flag) then\n"
-            "        b(widx1,widx2) = 0.0\n"
+            "  do widx2 = 1, SIZE(a(:,:), dim=2), 1\n"
+            "    do widx1 = 1, SIZE(a(:,:), dim=1), 1\n"
+            "      if (a(LBOUND(a, dim=1) + widx1 - 1,"
+            "LBOUND(a, dim=2) + widx2 - 1) < flag) then\n"
+            "        b(LBOUND(b, dim=1) + widx1 - 1,"
+            "LBOUND(b, dim=2) + widx2 - 1) = 0.0\n"
             "      end if\n"
             "    enddo\n"
             "  enddo\n"
@@ -281,10 +283,12 @@ def test_kernels_around_where_stmt(parser):
     new_code = str(psy.gen)
     assert ("  a(:,:) = 1.0\n"
             "  !$acc kernels\n"
-            "  do widx2 = 1, SIZE(a, 2), 1\n"
-            "    do widx1 = 1, SIZE(a, 1), 1\n"
-            "      if (a(widx1,widx2) < flag) then\n"
-            "        b(widx1,widx2) = 0.0\n"
+            "  do widx2 = 1, SIZE(a(:,:), dim=2), 1\n"
+            "    do widx1 = 1, SIZE(a(:,:), dim=1), 1\n"
+            "      if (a(LBOUND(a, dim=1) + widx1 - 1,"
+            "LBOUND(a, dim=2) + widx2 - 1) < flag) then\n"
+            "        b(LBOUND(b, dim=1) + widx1 - 1,"
+            "LBOUND(b, dim=2) + widx2 - 1) = 0.0\n"
             "      end if\n"
             "    enddo\n"
             "  enddo\n"
diff --git a/src/psyclone/tests/psyir/frontend/fparser2_test.py b/src/psyclone/tests/psyir/frontend/fparser2_test.py
index a54c3f668..c3a454314 100644
--- a/src/psyclone/tests/psyir/frontend/fparser2_test.py
+++ b/src/psyclone/tests/psyir/frontend/fparser2_test.py
@@ -1553,8 +1553,8 @@ def test_process_save_attribute_declarations(parser):
     assert isinstance(fake_parent.symbol_table.lookup("var4").interface,
                       StaticInterface)
 
-    # Test when is part of an UnknownDataType (target attribute in this case)
-    # it becomes an UnknownInterface
+    # Test that when it is part of an UnsupportedType (target attribute in
+    # this case) it becomes an UnknownInterface.
     reader = FortranStringReader("integer, target :: var5")
     fparser2spec = Specification_Part(reader).content[0]
     processor.process_declarations(fake_parent, [fparser2spec], [])
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 334b414b6..cbd36242b 100644
--- a/src/psyclone/tests/psyir/frontend/fparser2_where_handler_test.py
+++ b/src/psyclone/tests/psyir/frontend/fparser2_where_handler_test.py
@@ -45,9 +45,9 @@ from fparser.two import Fortran2003
 from psyclone.errors import InternalError
 from psyclone.psyir.frontend.fparser2 import Fparser2Reader
 from psyclone.psyir.nodes import (
-    Schedule, CodeBlock, Loop, ArrayReference, Assignment, Literal, Reference,
-    UnaryOperation, BinaryOperation, IfBlock, Call, Routine, Container, Range,
-    ArrayMember, StructureReference, IntrinsicCall)
+    ArrayMember, ArrayReference, Assignment, BinaryOperation, Call, CodeBlock,
+    Container, IfBlock, IntrinsicCall, Literal, Loop, Range, Routine, Schedule,
+    UnaryOperation)
 from psyclone.psyir.symbols import (
     DataSymbol, ArrayType, ScalarType, REAL_TYPE, INTEGER_TYPE,
     UnresolvedInterface)
@@ -327,7 +327,11 @@ def test_where_within_loop(fortran_reader):
     assert isinstance(where_loop.loop_body[0], IfBlock)
     assign = where_loop.loop_body[0].if_body[0]
     assert isinstance(assign, Assignment)
-    assert [idx.name for idx in assign.lhs.indices] == ["widx1", "jl"]
+    assert (assign.lhs.indices[0].debug_string() ==
+            "LBOUND(var2, dim=1) + widx1 - 1")
+    assert assign.lhs.indices[1].debug_string() == "jl"
+    assert where_loop.start_expr.value == "1"
+    assert where_loop.stop_expr.debug_string() == "SIZE(var, dim=1)"
 
 
 @pytest.mark.usefixtures("parser")
@@ -347,19 +351,16 @@ def test_basic_where():
         assert "was_where" in loop.annotations
         assert isinstance(loop.ast, Fortran2003.Where_Construct)
 
-    assert isinstance(loops[0].children[0], Literal)
-    assert isinstance(loops[0].children[1], IntrinsicCall)
-    assert loops[0].children[1].intrinsic == IntrinsicCall.Intrinsic.SIZE
-    assert str(loops[0].children[1].children[0]) == "Reference[name:'dry']"
+    assert isinstance(loops[0].start_expr, Literal)
+    assert loops[0].stop_expr.debug_string() == "SIZE(dry, dim=3)"
 
     ifblock = loops[2].loop_body[0]
     assert isinstance(ifblock, IfBlock)
     assert "was_where" in ifblock.annotations
-    assert ("ArrayReference[name:'dry']\n"
-            "Reference[name:'widx1']\n"
-            "Reference[name:'widx2']\n"
-            "Reference[name:'widx3']\n"
-            in str(ifblock.condition))
+    assert (ifblock.condition.debug_string() ==
+            "dry(LBOUND(dry, dim=1) + widx1 - 1,"
+            "LBOUND(dry, dim=2) + widx2 - 1,"
+            "LBOUND(dry, dim=3) + widx3 - 1)")
 
 
 @pytest.mark.usefixtures("parser")
@@ -384,9 +385,173 @@ def test_where_array_subsections():
     # Check that the array reference is indexed correctly
     assign = ifblock.if_body[0]
     assert isinstance(assign, Assignment)
-    assert isinstance(assign.lhs.children[0], Reference)
-    assert assign.lhs.children[0].name == "widx1"
-    assert assign.lhs.children[2].name == "widx2"
+    assert isinstance(assign.lhs.children[0], BinaryOperation)
+    assert (assign.lhs.children[0].debug_string() ==
+            "LBOUND(z1_st, dim=1) + widx1 - 1")
+    assert (assign.lhs.children[2].debug_string() ==
+            "LBOUND(z1_st, dim=3) + widx2 - 1")
+
+
+def test_where_mask_starting_value(fortran_reader, fortran_writer):
+    '''
+    Check handling of a case where the mask array is indexed from values other
+    than unity.
+
+    # TODO #949 - we can't currently take advantage of any knowledge of the
+    # declared lower bounds of arrays because the fparser2 frontend doesn't yet
+    # capture this information (we get an UnsupportedFortranType).
+    '''
+    code = '''\
+    program my_sub
+      use some_mod
+      real, dimension(-5:5,-5:5) :: picefr
+      real, DIMENSION(11,11,jpl) :: zevap_ice
+      real, dimension(-2:8,jpl,-3:7) :: snow
+      real, dimension(-22:0,jpl,-32:0) :: slush
+
+      WHERE( picefr(:,:) > 1.e-10 )
+        zevap_ice(:,:,1) = snow(:,3,:) * frcv(jpr_ievp)%z3(:,:,1) / picefr(:,:)
+      ELSEWHERE
+        zevap_ice(:,:,1) = snow(:,map(jpl),:) + slush(-22:-11,jpl,-32:-21)
+      END WHERE
+    end program my_sub
+'''
+    psyir = fortran_reader.psyir_from_source(code)
+    output = fortran_writer(psyir)
+    print(output)
+    expected = '''\
+  do widx2 = 1, SIZE(picefr, dim=2), 1
+    do widx1 = 1, SIZE(picefr, dim=1), 1
+      if (picefr(LBOUND(picefr, dim=1) + widx1 - 1,\
+LBOUND(picefr, dim=2) + widx2 - 1) > 1.e-10) then
+        zevap_ice(LBOUND(zevap_ice, dim=1) + widx1 - 1,\
+LBOUND(zevap_ice, dim=2) + widx2 - 1,1) = \
+snow(LBOUND(snow, dim=1) + widx1 - 1,3,LBOUND(snow, dim=3) + widx2 - 1) * \
+frcv(jpr_ievp)%z3(LBOUND(frcv(jpr_ievp)%z3, dim=1) + widx1 - 1,\
+LBOUND(frcv(jpr_ievp)%z3, dim=2) + widx2 - 1,1) / \
+picefr(LBOUND(picefr, dim=1) + widx1 - 1,LBOUND(picefr, dim=2) + widx2 - 1)
+      else
+        zevap_ice(LBOUND(zevap_ice, dim=1) + widx1 - 1,\
+LBOUND(zevap_ice, dim=2) + widx2 - 1,1) = \
+snow(LBOUND(snow, dim=1) + widx1 - 1,map(jpl),\
+LBOUND(snow, dim=3) + widx2 - 1) + slush(-22 + widx1 - 1,jpl,-32 + widx2 - 1)
+'''
+    assert expected in output
+
+
+def test_where_mask_is_slice(fortran_reader, fortran_writer):
+    '''
+    Check that the correct loop bounds and index expressions are created
+    when the mask expression uses a slice with specified bounds.
+    '''
+    code = '''\
+    program my_sub
+      use some_mod
+      WHERE( picefr(2:4,jstart:jstop) > 1.e-10 )
+        zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) / picefr(:,:)
+      ELSEWHERE ( picefr(1:3,jstart:jstop) > 4.e-10)
+        zevap_ice(:,:,1) = 0.0
+      END WHERE
+    end program my_sub
+'''
+    psyir = fortran_reader.psyir_from_source(code)
+    out = fortran_writer(psyir)
+    # Check that created loops have the correct number of iterations
+    assert "do widx2 = 1, jstop - jstart + 1, 1" in out
+    assert "do widx1 = 1, 4 - 2 + 1, 1" in out
+    # Check that the indexing into the mask expression uses the lower bounds
+    # specified in the original slice.
+    assert "if (picefr(2 + widx1 - 1,jstart + widx2 - 1) > 1.e-10)" in out
+    assert ("zevap_ice(LBOUND(zevap_ice, dim=1) + widx1 - 1,"
+            "LBOUND(zevap_ice, dim=2) + widx2 - 1,1) = 0.0" in out)
+    # If the lower bound of the slice is unity then we can use the loop
+    # index directly.
+    assert "if (picefr(widx1,jstart + widx2 - 1) > 4.e-10)" in out
+
+
+def test_where_body_containing_sum_with_dim(fortran_reader, fortran_writer):
+    '''
+    Since a SUM(x, dim=y) performs a reduction but produces an array, we need
+    to replace it with a temporary in order to translate the WHERE into
+    canonical form. We can't do that without being able to declare a suitable
+    temporary and that requires TODO #1799.
+    '''
+    code = '''\
+    module my_mod
+    contains
+    subroutine my_sub(picefr)
+      use some_mod
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   picefr
+      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zevap_ice
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple
+      WHERE( picefr(:,:) > 1.e-10 )
+        zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) * &
+                           SUM( a_i_last_couple, dim=3 ) / picefr(:,:)
+      ELSEWHERE
+        zevap_ice(:,:,1) = 0.0
+      END WHERE
+    end subroutine my_sub
+    end module my_mod
+'''
+    psyir = fortran_reader.psyir_from_source(code)
+    routine = psyir.walk(Routine)[0]
+    assert isinstance(routine[0], CodeBlock)
+    output = fortran_writer(psyir)
+    assert "WHERE (picefr" in output
+    assert "SUM(a_i_last_couple, dim = 3) / picefr(:, :)" in output
+
+
+def test_where_containing_sum_no_dim(fortran_reader, fortran_writer):
+    '''
+    Since a SUM without a dim argument always produces a scalar we can
+    translate a WHERE containing it into canonical form.
+    '''
+    code = '''\
+    module my_mod
+    contains
+    subroutine my_sub(picefr)
+      use some_mod
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   picefr
+      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zevap_ice
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple
+      WHERE( picefr(:,:) > SUM(picefr) )
+        zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) * &
+                           SUM( a_i_last_couple ) / picefr(:,:)
+      ELSEWHERE
+        zevap_ice(:,:,1) = 0.0
+      END WHERE
+    end subroutine my_sub
+    end module my_mod
+'''
+    psyir = fortran_reader.psyir_from_source(code)
+    routine = psyir.walk(Routine)[0]
+    assert isinstance(routine[0], Loop)
+    output = fortran_writer(psyir)
+    assert ("SUM(a_i_last_couple) / picefr(LBOUND(picefr, dim=1) + widx1 - 1,"
+            "LBOUND(picefr, dim=2) + widx2 - 1)" in output)
+
+
+def test_where_mask_containing_sum_with_dim(fortran_reader):
+    '''Since a SUM(x, dim=y) appearing in a mask expression performs
+    a reduction but produces an array, we need to replace it with a
+    temporary in order to translate the WHERE into canonical form
+    (since the mask expression determines the number of nested loops
+    required). We can't do that without being able to declare a
+    suitable temporary and that requires TODO #1799.
+
+    '''
+    code = '''\
+    subroutine my_sub(v2)
+      REAL(wp), INTENT(in), DIMENSION(:,:) :: v2
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: v3
+      where(sum(v2(:,:), dim=2) > 0.0)
+        v3(:) = 1.0
+      end where
+    end subroutine my_sub
+'''
+    psyir = fortran_reader.psyir_from_source(code)
+    routine = psyir.walk(Routine)[0]
+    assert isinstance(routine[0], CodeBlock)
 
 
 @pytest.mark.usefixtures("parser")
@@ -453,9 +618,10 @@ def test_elsewhere():
     assert isinstance(ifblock.ast, Fortran2003.Where_Construct)
     assert isinstance(ifblock.condition, BinaryOperation)
     assert ifblock.condition.operator == BinaryOperation.Operator.GT
-    assert ("ArrayReference[name:'ptsu']\n"
-            "Reference[name:'widx1']\n" in str(ifblock.condition.children[0]))
-    assert "Literal[value:'10." in str(ifblock.condition.children[1])
+    assert (ifblock.condition.debug_string() ==
+            "ptsu(LBOUND(ptsu, dim=1) + widx1 - 1,"
+            "LBOUND(ptsu, dim=2) + widx2 - 1,"
+            "LBOUND(ptsu, dim=3) + widx3 - 1) > 10._wp")
     # Check that this IF block has an else body which contains another IF
     assert ifblock.else_body is not None
     ifblock2 = ifblock.else_body[0]
@@ -463,8 +629,10 @@ def test_elsewhere():
     assert isinstance(ifblock2, IfBlock)
     assert isinstance(ifblock2.condition, BinaryOperation)
     assert ifblock2.condition.operator == BinaryOperation.Operator.LT
-    assert ("ArrayReference[name:'ptsu']\n"
-            "Reference[name:'widx1']\n" in str(ifblock2.condition.children[0]))
+    assert (ifblock2.condition.debug_string() ==
+            "ptsu(LBOUND(ptsu, dim=1) + widx1 - 1,"
+            "LBOUND(ptsu, dim=2) + widx2 - 1,"
+            "LBOUND(ptsu, dim=3) + widx3 - 1) < 0.0_wp")
     # Check that this IF block too has an else body
     assert isinstance(ifblock2.else_body[0], Assignment)
     # Check that we have three assignments of the correct form and with the
@@ -473,12 +641,10 @@ def test_elsewhere():
     assert len(assigns) == 3
     for assign in assigns:
         assert isinstance(assign.lhs, ArrayReference)
-        refs = assign.lhs.walk(Reference)
-        assert len(refs) == 4
-        assert refs[1:4] == assign.lhs.indices
-        assert assign.lhs.name == "z1_st"
-        assert ([idx.name for idx in assign.lhs.indices] ==
-                ["widx1", "widx2", "widx3"])
+        assert (assign.lhs.debug_string() ==
+                "z1_st(LBOUND(z1_st, dim=1) + widx1 - 1,"
+                "LBOUND(z1_st, dim=2) + widx2 - 1,"
+                "LBOUND(z1_st, dim=3) + widx3 - 1)")
         assert isinstance(assign.parent.parent, IfBlock)
 
     assert isinstance(assigns[0].rhs, BinaryOperation)
@@ -526,6 +692,26 @@ def test_where_stmt():
     assert isinstance(fake_parent[0], Loop)
 
 
+def test_where_stmt_no_reduction(fortran_reader, fortran_writer):
+    '''
+    Test that a WHERE statement containing an intrinsic reduction is put
+    into a CodeBlock.
+    '''
+    code = '''\
+    program my_prog
+    integer, dimension(10,10,5) :: a3d
+    integer, dimension(10,10) :: a_i, at_i, rn_amax_2d
+    WHERE( at_i(:,:) > rn_amax_2d(:,:) ) a_i(:,:,jl) = a_i(:,:,jl) * \
+rn_amax_2d(:,:) / sum(a3d, dim=3)
+    end program my_prog'''
+    psyir = fortran_reader.psyir_from_source(code)
+    routine = psyir.walk(Routine)[0]
+    assert isinstance(routine[0], CodeBlock)
+    output = fortran_writer(psyir)
+    assert ("WHERE (at_i(:, :) > rn_amax_2d(:, :)) a_i(:, :, jl) = "
+            "a_i(:, :, jl) * rn_amax_2d(:, :) / SUM(a3d, dim = 3)" in output)
+
+
 def test_where_ordering(parser):
     ''' Check that the generated schedule has the correct ordering when
     a WHERE construct is processed.
@@ -588,13 +774,11 @@ def test_where_derived_type(fortran_reader, fortran_writer, code, size_arg):
     loops = psyir.walk(Loop)
     assert len(loops) == 2
     assert isinstance(loops[1].stop_expr, IntrinsicCall)
-    assert loops[1].stop_expr.intrinsic == IntrinsicCall.Intrinsic.SIZE
-    assert isinstance(loops[1].stop_expr.children[0], StructureReference)
-    assert fortran_writer(loops[1].stop_expr.children[0]) == size_arg
+    assert loops[1].stop_expr.debug_string() == f"SIZE({size_arg}, dim=1)"
     assert isinstance(loops[1].loop_body[0], IfBlock)
     # All Range nodes should have been replaced
     assert not loops[0].walk(Range)
     # All ArrayMember accesses should now use the `widx1` loop variable
     array_members = loops[0].walk(ArrayMember)
     for member in array_members:
-        assert member.indices[0].symbol.name == "widx1"
+        assert "+ widx1 - 1" in member.indices[0].debug_string()
diff --git a/src/psyclone/tests/psyir/nodes/array_mixin_test.py b/src/psyclone/tests/psyir/nodes/array_mixin_test.py
index afade7093..632fc991d 100644
--- a/src/psyclone/tests/psyir/nodes/array_mixin_test.py
+++ b/src/psyclone/tests/psyir/nodes/array_mixin_test.py
@@ -488,10 +488,22 @@ def test_aref_get_full_range_unknown_size(extent):
             "less than the number of dimensions '2'." in str(excinfo.value))
 
 
+# _extent
+
+def test_arraymixin_extent(fortran_reader):
+    '''Tests for the _extent() method.'''
+    atype = ArrayType(INTEGER_TYPE, [10])
+    asym = DataSymbol("array", atype)
+    # Access to single array element just has extent of 1.
+    ref = ArrayReference.create(asym, [Literal("1", INTEGER_TYPE)])
+    assert ref._extent(0).value == "1"
+    ref2 = ArrayReference.create(asym, [ref.get_full_range(0)])
+    assert ref2._extent(0).debug_string() == "10"
+
 # _get_effective_shape
 
 
-def test_get_effective_shape(fortran_reader, fortran_writer):
+def test_get_effective_shape(fortran_reader):
     '''Tests for the _get_effective_shape() method.'''
     code = (
         "subroutine test()\n"
@@ -499,6 +511,12 @@ def test_get_effective_shape(fortran_reader, fortran_writer):
         "  integer :: indices(8,3)\n"
         "  real a(10), b(10,10)\n"
         "  a(1:10) = 0.0\n"
+        "  a(1:10:3) = 0.0\n"
+        "  a(:) = 0.0\n"
+        "  a(2:) = 0.0\n"
+        "  a(:5) = 0.0\n"
+        "  a(::4) = 0.0\n"
+        "  a(lbound(b,1):ubound(b,2)) = 0.0\n"
         "  b(indices(2:3,1), 2:5) = 2.0\n"
         "  b(indices(2:3,1:2), 2:5) = 2.0\n"
         "  a(f()) = 2.0\n"
@@ -507,33 +525,71 @@ def test_get_effective_shape(fortran_reader, fortran_writer):
     psyir = fortran_reader.psyir_from_source(code)
     routine = psyir.walk(Routine)[0]
     # Direct array slice.
-    shape = routine.children[0].lhs._get_effective_shape()
+    child_idx = 0
+    shape = routine.children[child_idx].lhs._get_effective_shape()
     assert len(shape) == 1
-    assert isinstance(shape[0], BinaryOperation)
+    assert isinstance(shape[0], Literal)
+    assert shape[0].value == "10"
+    # Array slice with non-unit step.
+    child_idx += 1
+    shape = routine.children[child_idx].lhs._get_effective_shape()
+    assert len(shape) == 1
+    assert shape[0].debug_string() == "(10 - 1) / 3 + 1"
+    # Full array slice without bounds.
+    child_idx += 1
+    shape = routine.children[child_idx].lhs._get_effective_shape()
+    assert len(shape) == 1
+    assert "SIZE(a, dim=1)" in shape[0].debug_string()
+    # Array slice with only lower-bound specified.
+    child_idx += 1
+    shape = routine.children[child_idx].lhs._get_effective_shape()
+    assert len(shape) == 1
+    assert shape[0].debug_string() == "UBOUND(a, dim=1) - 2 + 1"
+    # Array slice with only upper-bound specified.
+    child_idx += 1
+    shape = routine.children[child_idx].lhs._get_effective_shape()
+    assert len(shape) == 1
+    assert shape[0].debug_string() == "5 - LBOUND(a, dim=1) + 1"
+    # Array slice with only step specified.
+    child_idx += 1
+    shape = routine.children[child_idx].lhs._get_effective_shape()
+    assert len(shape) == 1
+    # Since step is not unity, this is not a 'full-range' access so we still
+    # end up with UBOUND and LBOUND, even though SIZE would be nicer.
+    assert (shape[0].debug_string() ==
+            "(UBOUND(a, dim=1) - LBOUND(a, dim=1)) / 4 + 1")
+    # Array slice defined using LBOUND and UBOUND intrinsics but for a
+    # different array altogether.
+    child_idx += 1
+    shape = routine.children[child_idx].lhs._get_effective_shape()
+    assert len(shape) == 1
+    assert shape[0].debug_string() == "UBOUND(b, 2) - LBOUND(b, 1) + 1"
     # Indirect array slice.
-    shape = routine.children[1].lhs._get_effective_shape()
+    child_idx += 1
+    shape = routine.children[child_idx].lhs._get_effective_shape()
     assert len(shape) == 2
     assert isinstance(shape[0], BinaryOperation)
-    code = fortran_writer(shape[0])
     # An ArrayType does not store the number of elements, just lower and upper
     # bounds. Therefore, we end up recursively computing the no. of elements.
     # The answer is still "2" though!
-    assert code == "((3 - 2) / 1 + 1 - 1) / 1 + 1"
-    code = fortran_writer(shape[1])
-    assert code == "(5 - 2) / 1 + 1"
+    assert shape[0].debug_string() == "3 - 2 + 1"
+    assert shape[1].debug_string() == "5 - 2 + 1"
     # An indirect array slice can only be 1D.
+    child_idx += 1
     with pytest.raises(InternalError) as err:
-        _ = routine.children[2].lhs._get_effective_shape()
+        _ = routine.children[child_idx].lhs._get_effective_shape()
     assert ("array defining a slice of a dimension of another array must be "
             "1D but 'indices' used to index into 'b' has 2 dimensions" in
             str(err.value))
     # Indirect array access using function call.
+    child_idx += 1
     with pytest.raises(NotImplementedError) as err:
-        _ = routine.children[3].lhs._get_effective_shape()
+        _ = routine.children[child_idx].lhs._get_effective_shape()
     assert "include a function call or expression" in str(err.value)
     # Array access with expression in indices.
+    child_idx += 1
     with pytest.raises(NotImplementedError) as err:
-        _ = routine.children[4].lhs._get_effective_shape()
+        _ = routine.children[child_idx].lhs._get_effective_shape()
     assert "include a function call or expression" in str(err.value)
 
 
diff --git a/src/psyclone/tests/psyir/nodes/array_reference_test.py b/src/psyclone/tests/psyir/nodes/array_reference_test.py
index 40a1e8d71..cfd72f90e 100644
--- a/src/psyclone/tests/psyir/nodes/array_reference_test.py
+++ b/src/psyclone/tests/psyir/nodes/array_reference_test.py
@@ -43,11 +43,12 @@ import pytest
 from psyclone.errors import GenerationError, InternalError
 from psyclone.psyir.backend.fortran import FortranWriter
 from psyclone.psyir.nodes.node import colored
-from psyclone.psyir.nodes import Reference, ArrayReference, Assignment, \
-    Literal, BinaryOperation, Range, KernelSchedule, IntrinsicCall
+from psyclone.psyir.nodes import (
+    Reference, ArrayReference, Assignment,
+    Literal, BinaryOperation, Range, KernelSchedule, IntrinsicCall)
 from psyclone.psyir.symbols import (
     ArrayType, DataSymbol, DataTypeSymbol, UnresolvedType, ScalarType,
-    REAL_SINGLE_TYPE, INTEGER_SINGLE_TYPE, REAL_TYPE, INTEGER_TYPE,
+    REAL_SINGLE_TYPE, INTEGER_SINGLE_TYPE, REAL_TYPE, Symbol, INTEGER_TYPE,
     UnsupportedFortranType)
 from psyclone.tests.utilities import check_links
 
@@ -506,6 +507,7 @@ def test_array_datatype():
     one = Literal("1", INTEGER_TYPE)
     two = Literal("2", INTEGER_TYPE)
     four = Literal("4", INTEGER_TYPE)
+    six = Literal("6", INTEGER_TYPE)
     # Reference to a single element of an array.
     aref = ArrayReference.create(test_sym, [one])
     assert aref.datatype == REAL_TYPE
@@ -523,7 +525,20 @@ def test_array_datatype():
     assert isinstance(upper, BinaryOperation)
     # The easiest way to check the expression is to use debug_string()
     code = upper.debug_string()
-    assert code == "(4 - 2) / 1 + 1"
+    assert code == "4 - 2 + 1"
+    # Reference to a non-contiguous 1D sub-array of a 2D array.
+    ucref = ArrayReference.create(test_sym2d, [two.copy(),
+                                               Range.create(two.copy(),
+                                                            six.copy(),
+                                                            two.copy())])
+    assert isinstance(ucref.datatype, ArrayType)
+    assert ucref.datatype.intrinsic == ScalarType.Intrinsic.REAL
+    assert len(ucref.datatype.shape) == 1
+    # The sub-array will have a lower bound of one.
+    assert ucref.datatype.shape[0].lower == one
+    # Upper bound must be computed as (stop - start)/step + 1
+    upper = ucref.datatype.shape[0].upper
+    assert upper.debug_string() == "(6 - 2) / 2 + 1"
     # Reference to a single element of an array of structures.
     stype = DataTypeSymbol("grid_type", UnresolvedType())
     atype = ArrayType(stype, [10])
@@ -531,23 +546,23 @@ def test_array_datatype():
     aref = ArrayReference.create(asym, [two.copy()])
     assert aref.datatype is stype
     # Reference to a single element of an array of UnsupportedType.
-    unsuppored_sym = DataSymbol(
-        "unsuppored",
-        UnsupportedFortranType("real, dimension(5), pointer :: unsuppored"))
-    aref = ArrayReference.create(unsuppored_sym, [two.copy()])
+    unsupported_sym = DataSymbol(
+        "unsupported",
+        UnsupportedFortranType("real, dimension(5), pointer :: unsupported"))
+    aref = ArrayReference.create(unsupported_sym, [two.copy()])
     assert isinstance(aref.datatype, UnresolvedType)
     # Reference to a single element of an array of UnsupportedType but with
     # partial type information.
-    not_quite_unsuppored_sym = DataSymbol(
-        "unsuppored",
+    not_quite_unsupported_sym = DataSymbol(
+        "unsupported",
         UnsupportedFortranType(
-            "real, dimension(5), pointer :: unsuppored",
+            "real, dimension(5), pointer :: unsupported",
             partial_datatype=ArrayType(REAL_SINGLE_TYPE, [5])))
-    bref = ArrayReference.create(not_quite_unsuppored_sym, [two.copy()])
+    bref = ArrayReference.create(not_quite_unsupported_sym, [two.copy()])
     assert bref.datatype == REAL_SINGLE_TYPE
     # A sub-array of UnsupportedFortranType.
     aref3 = ArrayReference.create(
-                unsuppored_sym, [Range.create(two.copy(), four.copy())])
+                unsupported_sym, [Range.create(two.copy(), four.copy())])
     # We know the result is an ArrayType
     assert isinstance(aref3.datatype, ArrayType)
     assert aref3.datatype.shape[0].lower == one
@@ -557,8 +572,20 @@ def test_array_datatype():
     assert isinstance(aref3.datatype.intrinsic, UnresolvedType)
     # A whole array of UnsupportedType should simply have the same datatype as
     # the original symbol.
-    aref4 = ArrayReference.create(not_quite_unsuppored_sym, [":"])
-    assert aref4.datatype == not_quite_unsuppored_sym.datatype
+    aref4 = ArrayReference.create(not_quite_unsupported_sym, [":"])
+    assert aref4.datatype == not_quite_unsupported_sym.datatype
+    # When the Reference is just to a Symbol. Although `create` forbids this,
+    # it is possible for the fparser2 frontend to create such a construct.
+    generic_sym = Symbol("existential")
+    aref5 = ArrayReference(generic_sym)
+    aref5.addchild(two.copy())
+    assert isinstance(aref5.datatype, UnresolvedType)
+    aref5.addchild(Range.create(two.copy(), four.copy()))
+    dtype5 = aref5.datatype
+    assert isinstance(dtype5, ArrayType)
+    assert isinstance(dtype5.intrinsic, UnresolvedType)
+    assert dtype5.shape[0].lower.value == "1"
+    assert dtype5.shape[0].upper.debug_string() == "4 - 2 + 1"
 
 
 def test_array_create_colon(fortran_writer):
diff --git a/src/psyclone/tests/psyir/nodes/operation_test.py b/src/psyclone/tests/psyir/nodes/operation_test.py
index 76b692278..3f3d169ea 100644
--- a/src/psyclone/tests/psyir/nodes/operation_test.py
+++ b/src/psyclone/tests/psyir/nodes/operation_test.py
@@ -331,7 +331,7 @@ def test_binaryop_array_section_datatype():
     dtype1 = binaryoperation.datatype
     assert len(dtype1.shape) == 1
     assert dtype1.shape[0].lower.value == "1"
-    assert dtype1.shape[0].upper.debug_string() == "(3 - 2) / 1 + 1"
+    assert dtype1.shape[0].upper.debug_string() == "3 - 2 + 1"
     # Repeat but for a non-contiguous 1D section of a rank 2, integer array.
     ref3 = ArrayReference.create(
         DataSymbol("tmp3", iarrtype),
@@ -343,7 +343,7 @@ def test_binaryop_array_section_datatype():
     assert dtype3.intrinsic == REAL_SINGLE_TYPE.intrinsic
     assert len(dtype3.shape) == 1
     assert dtype3.shape[0].lower.value == "1"
-    assert dtype3.shape[0].upper.debug_string() == "(3 - 2) / 1 + 1"
+    assert dtype3.shape[0].upper.debug_string() == "3 - 2 + 1"
 
 
 def test_binaryop_datatype_recursion():
@@ -371,7 +371,7 @@ def test_binaryop_datatype_recursion():
     assert dtype1.intrinsic == REAL_SINGLE_TYPE.intrinsic
     assert len(dtype1.shape) == 1
     assert dtype1.shape[0].lower.value == "1"
-    assert dtype1.shape[0].upper.debug_string() == "(3 - 2) / 1 + 1"
+    assert dtype1.shape[0].upper.debug_string() == "3 - 2 + 1"
 
 
 def test_binaryop_structure_datatype():