Repository: stfc/PSyclone

Issue:
Fortran WHERE construct unsupported
Hello,
processing some legacy FORTRAN (UKCA glomap) that has several (many) WHERE statements; there are instances which PSyclone cannot support. i.e. when an array (1:nbox) is coded without the array index subscript... 
The example shown is declared as corrh(1:nbox), is occasionally written as corrh. I wrote a toy sample to demo the situation.
Also noted that the write_stmnt comment is interleaved with the WHERE unsupported report. Line 27 is the actual troublesome statement.
[corrh_where_f90.txt](https://github.com/user-attachments/files/16494639/corrh_where_f90.txt)

==== quoted  ===
21   ! PSyclone CodeBlock (unsupported code) reason:
 22   !  - Unsupported statement: Write_Stmt
 23   !  - Only WHERE constructs using explicit array notation (e.g. my_array(:,    :)) are supported but found 
 24 !& '[Level_4_Expr(Name('corrh'), '<', Real_Literal_Constant('1.0', None))]'.
 25   !  - Unsupported statement: Write_Stmt
 26   WRITE(6, '(20F5.2,2x)') corrh(1 : 20)
 27   WHERE (corrh < 1.0) corrh = 1.05

=== end quoted ===

Discussion:

Patch:
diff --git a/changelog b/changelog
index f730cbc74..a8f79b45f 100644
--- a/changelog
+++ b/changelog
@@ -228,6 +228,9 @@
 	78) PR #2678 for #2636. Alter LFRic PSy-layer to lookup nlevels for
 	each kernel
 
+	79) PR #2687 for #2684. Add support for WHERE constructs that contain
+	references without range-notation (and prevent some incorrect cases).
+
 release 2.5.0 14th of February 2024
 
 	1) PR #2199 for #2189. Fix bugs with missing maps in enter data
diff --git a/examples/nemo/scripts/acc_kernels_trans.py b/examples/nemo/scripts/acc_kernels_trans.py
index 3099bc87a..a08439130 100755
--- a/examples/nemo/scripts/acc_kernels_trans.py
+++ b/examples/nemo/scripts/acc_kernels_trans.py
@@ -63,7 +63,7 @@ from psyclone.errors import InternalError
 from psyclone.psyGen import TransInfo
 from psyclone.psyir.nodes import (
     IfBlock, ArrayReference, Assignment, BinaryOperation, Loop, Routine,
-    Literal, Call, ACCLoopDirective)
+    Literal, ACCLoopDirective)
 from psyclone.psyir.transformations import (ACCKernelsTrans, ACCUpdateTrans,
                                             TransformationError, ProfileTrans)
 from psyclone.transformations import ACCEnterDataTrans
@@ -392,17 +392,12 @@ def trans(psyir):
     for subroutine in psyir.walk(Routine):
         print(f"Transforming subroutine: {subroutine.name}")
 
-        # In the lib_fortran file we annotate each routine that does not
-        # have a Loop or unsupported Calls with the OpenACC Routine Directive
+        # In the lib_fortran file we annotate each routine of the SIGN_*
+        # interface with the OpenACC Routine Directive
         if psyir.name == "lib_fortran.f90":
-            if not subroutine.walk(Loop):
-                calls = subroutine.walk(Call)
-                if all(call.is_available_on_device() for call in calls):
-                    # SIGN_ARRAY_1D has a CodeBlock because of a WHERE without
-                    # array notation. (TODO #717)
-                    ACC_ROUTINE_TRANS.apply(subroutine,
-                                            options={"force": True})
-                    continue
+            if subroutine.name.lower().startswith("sign_"):
+                ACC_ROUTINE_TRANS.apply(subroutine)
+                continue
 
         # Attempt to add OpenACC directives unless we are ignoring this routine
         if subroutine.name.lower() not in ACC_IGNORE:
diff --git a/examples/nemo/scripts/acc_loops_trans.py b/examples/nemo/scripts/acc_loops_trans.py
index bf7bb9994..6d293d4b6 100755
--- a/examples/nemo/scripts/acc_loops_trans.py
+++ b/examples/nemo/scripts/acc_loops_trans.py
@@ -40,9 +40,9 @@ directives into Nemo code. '''
 from utils import (
     insert_explicit_loop_parallelism, normalise_loops, add_profiling,
     enhance_tree_information, NOT_PERFORMANT, NOT_WORKING)
-from psyclone.psyir.nodes import Loop, Routine
+from psyclone.psyir.nodes import Routine
 from psyclone.transformations import (
-    ACCParallelTrans, ACCLoopTrans, ACCRoutineTrans, TransformationError)
+    ACCParallelTrans, ACCLoopTrans, ACCRoutineTrans)
 
 PROFILING_ENABLED = True
 
@@ -108,17 +108,12 @@ def trans(psyir):
                 hoist_expressions=True
         )
 
-        # For performance in lib_fortran, mark serial routines as GPU-enabled
+        # In the lib_fortran file we annotate each routine of the SIGN_*
+        # interface with the OpenACC Routine Directive
         if psyir.name == "lib_fortran.f90":
-            if not subroutine.walk(Loop):
-                try:
-                    # We need the 'force' option.
-                    # SIGN_ARRAY_1D has a CodeBlock because of a WHERE without
-                    # array notation. (TODO #717)
-                    ACCRoutineTrans().apply(subroutine,
-                                            options={"force": True})
-                except TransformationError as err:
-                    print(err)
+            if subroutine.name.lower().startswith("sign_"):
+                ACCRoutineTrans().apply(subroutine)
+                continue
 
         insert_explicit_loop_parallelism(
             subroutine,
diff --git a/examples/nemo/scripts/omp_gpu_trans.py b/examples/nemo/scripts/omp_gpu_trans.py
index c890c8fe3..916dc9868 100755
--- a/examples/nemo/scripts/omp_gpu_trans.py
+++ b/examples/nemo/scripts/omp_gpu_trans.py
@@ -98,17 +98,14 @@ def trans(psyir):
                 hoist_expressions=True
         )
 
-        # For performance in lib_fortran, mark serial routines as GPU-enabled
+        # In the lib_fortran file we annotate each routine of the SIGN_*
+        # interface with the OpenMP Declare Target Directive
         if psyir.name == "lib_fortran.f90":
-            if not subroutine.walk(Loop):
-                try:
-                    # We need the 'force' option.
-                    # SIGN_ARRAY_1D has a CodeBlock because of a WHERE without
-                    # array notation. (TODO #717)
-                    OMPDeclareTargetTrans().apply(subroutine,
-                                                  options={"force": True})
-                except TransformationError as err:
-                    print(err)
+            if subroutine.name.lower().startswith("sign_"):
+                OMPDeclareTargetTrans().apply(subroutine)
+                # We continue parallelising inside the routine, but this could
+                # change if the parallelisation directive are not nestable, in
+                # which case we could add a 'continue' here
 
         # For now this is a special case for stpctl.f90 because it forces
         # loops to parallelise without many safety checks
diff --git a/src/psyclone/psyir/frontend/fparser2.py b/src/psyclone/psyir/frontend/fparser2.py
index 1bcfca9c7..0fdcb32a6 100644
--- a/src/psyclone/psyir/frontend/fparser2.py
+++ b/src/psyclone/psyir/frontend/fparser2.py
@@ -60,8 +60,6 @@ from psyclone.psyir.nodes import (
     Reference, Return, Routine, Schedule, StructureReference, UnaryOperation,
     WhileLoop)
 from psyclone.psyir.nodes.array_mixin import ArrayMixin
-from psyclone.psyir.nodes.array_of_structures_mixin import (
-    ArrayOfStructuresMixin)
 from psyclone.psyir.symbols import (
     ArgumentInterface, ArrayType, AutomaticInterface, CHARACTER_TYPE,
     CommonBlockInterface, ContainerSymbol, DataSymbol, DataTypeSymbol,
@@ -4179,109 +4177,48 @@ class Fparser2Reader():
                 parent.addchild(call)
                 call.children.extend(fake_parent.pop_all_children())
 
-    @staticmethod
-    def _array_notation_rank(node):
-        '''Check that the supplied candidate array reference uses supported
-        array notation syntax and return the rank of the sub-section
-        of the array that uses array notation. e.g. for a reference
-        "a(:, 2, :)" the rank of the sub-section is 2.
-
-        :param node: the reference to check.
-        :type node: :py:class:`psyclone.psyir.nodes.ArrayReference` or \
-            :py:class:`psyclone.psyir.nodes.ArrayMember` or \
-            :py:class:`psyclone.psyir.nodes.StructureReference`
-
-        :returns: rank of the sub-section of the array.
-        :rtype: int
-
-        :raises InternalError: if no ArrayMixin node with at least one \
-                               Range in its indices is found.
-        :raises InternalError: if two or more part references in a \
-                               structure reference contain ranges.
-        :raises NotImplementedError: if the supplied node is not of a \
-                                     supported type.
-        :raises NotImplementedError: if any ranges are encountered that are \
-                                     not for the full extent of the dimension.
-        '''
-        if isinstance(node, (ArrayReference, ArrayMember)):
-            array = node
-        elif isinstance(node, StructureReference):
-            array = None
-            arrays = node.walk((ArrayMember, ArrayOfStructuresMixin))
-            for part_ref in arrays:
-                if any(isinstance(idx, Range) for idx in part_ref.indices):
-                    if array:
-                        # Cannot have two or more part references that contain
-                        # ranges - this is not valid Fortran.
-                        raise InternalError(
-                            f"Found a structure reference containing two or "
-                            f"more part references that have ranges: "
-                            f"'{node.debug_string()}'. This is not valid "
-                            f"within a WHERE in Fortran.")
-                    array = part_ref
-            if not array:
-                raise InternalError(
-                    f"No array access found in node '{node.name}'")
-        else:
-            # This will result in a CodeBlock.
-            raise NotImplementedError(
-                f"Expected either an ArrayReference, ArrayMember or a "
-                f"StructureReference but got '{type(node).__name__}'")
-
-        # Only array refs using basic colon syntax are currently
-        # supported e.g. (a(:,:)).  Each colon is represented in the
-        # PSyIR as a Range node with first argument being an lbound
-        # binary operator, the second argument being a ubound operator
-        # and the third argument being an integer Literal node with
-        # value 1 i.e. a(:,:) is represented as
-        # a(lbound(a,1):ubound(a,1):1,lbound(a,2):ubound(a,2):1) in
-        # the PSyIR.
-        num_colons = 0
-        for idx_node in array.indices:
-            if isinstance(idx_node, Range):
-                # Found array syntax notation. Check that it is the
-                # simple ":" format.
-                if not array.is_full_range(array.index_of(idx_node)):
-                    raise NotImplementedError(
-                        "Only array notation of the form my_array(:, :, ...) "
-                        "is supported.")
-                num_colons += 1
-        return num_colons
-
     def _array_syntax_to_indexed(self, parent, loop_vars):
         '''
-        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. This indexing
-        is always done relative to the declared lower bound of the array being
-        accessed.
+        Utility function that modifies bare References to arrays as
+        ArrayReferences, and each ArrayReference that is not inside an
+        elemental call to use the provided loop index.
 
         :param parent: root of PSyIR sub-tree to search for Array
                        references to modify.
-        :type parent:  :py:class:`psyclone.psyir.nodes.Node`
+        :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
                                      found.
         '''
-        assigns = parent.walk(Assignment)
-        # Check that the LHS of any assignment uses array notation.
-        # Note that this will prevent any WHERE blocks that contain scalar
-        # assignments from being handled but is a necessary limitation until
-        # #717 is done and we interrogate the type of each symbol.
-        for assign in assigns:
-            _ = self._array_notation_rank(assign.lhs)
-
-        # TODO #717 if the supplied code accidentally omits array
-        # notation for an array reference on the RHS then we will
-        # identify it as a scalar and the code produced from the
-        # 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.
+        # pylint: disable=import-outside-toplevel
+        from psyclone.psyir.transformations import (
+                Reference2ArrayRangeTrans, TransformationError)
+        # Convert References to arrays to use the array range notation unless
+        # they have an IntrinsicCall parent.
+        for ref in parent.walk(Reference):
+            if isinstance(ref.symbol.interface, (ImportInterface,
+                                                 UnresolvedInterface)):
+                raise NotImplementedError(
+                        "PSyclone doesn't yet support reference to imported "
+                        "symbols inside WHERE clauses.")
+            call_ancestor = ref.ancestor(Call)
+            if (isinstance(ref.symbol, DataSymbol) and not call_ancestor):
+                try:
+                    Reference2ArrayRangeTrans().apply(ref)
+                except TransformationError:
+                    pass
+            elif (isinstance(ref.symbol, DataSymbol) and call_ancestor
+                  is not None and call_ancestor.is_elemental):
+                try:
+                    Reference2ArrayRangeTrans().apply(ref)
+                except TransformationError:
+                    pass
         table = parent.scope.symbol_table
         one = Literal("1", INTEGER_TYPE)
         arrays = parent.walk(ArrayMixin)
+
         first_rank = None
         for array in arrays:
             # Check that this is a supported array reference and that
@@ -4349,7 +4286,11 @@ class Fparser2Reader():
                     expr = BinaryOperation.create(
                         add_op, lbound, Reference(symbol))
                     expr2 = BinaryOperation.create(sub_op, expr, one.copy())
-                array.children[idx] = expr2
+                # Indices of an AoSReference start at index 1
+                if isinstance(array, ArrayOfStructuresReference):
+                    array.children[idx+1] = expr2
+                else:
+                    array.children[idx] = expr2
                 range_idx += 1
 
     def _where_construct_handler(self, node, parent):
@@ -4386,6 +4327,10 @@ class Fparser2Reader():
             not use array notation.
 
         '''
+        # pylint: disable=import-outside-toplevel
+        from psyclone.psyir.transformations import (
+                Reference2ArrayRangeTrans, TransformationError)
+
         def _contains_intrinsic_reduction(pnodes):
             '''
             Utility to check for Fortran intrinsics that perform a reduction
@@ -4474,17 +4419,25 @@ class Fparser2Reader():
         # parent for this logical expression we will repeat the processing.
         fake_parent = Assignment(parent=parent)
         self.process_nodes(fake_parent, logical_expr)
-        arrays = fake_parent.walk(ArrayMixin)
+        # We want to convert all the plain references that are arrays to use
+        # explicit array syntax.
+        # TODO 2722: Should have the same logic as array_syntax_to_indexed
+        # regarding UnresolvedInterface and Elemental calls?
+        references = fake_parent.walk(Reference)
+        for ref in references:
+            if isinstance(ref.symbol.interface, ImportInterface):
+                raise NotImplementedError(
+                        "PSyclone doesn't yet support reference to imported "
+                        "symbols inside WHERE clauses.")
+            intrinsic_ancestor = ref.ancestor(IntrinsicCall)
+            if (isinstance(ref.symbol, DataSymbol) and
+                    not intrinsic_ancestor):
+                try:
+                    Reference2ArrayRangeTrans().apply(ref)
+                except TransformationError:
+                    pass
 
-        if not arrays:
-            # If the PSyIR doesn't contain any Arrays then that must be
-            # 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 #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}'.")
+        arrays = fake_parent.walk(ArrayMixin)
 
         for array in arrays:
             if any(isinstance(idx, Range) for idx in array.indices):
@@ -4568,7 +4521,6 @@ 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)
@@ -4591,7 +4543,6 @@ class Fparser2Reader():
         # Now construct the body of the IF using the body of the WHERE
         sched = Schedule(parent=ifblock)
         ifblock.addchild(sched)
-
         if was_single_stmt:
             # We only had a single-statement WHERE
             self.process_nodes(sched, node.items[1:])
@@ -4663,7 +4614,6 @@ class Fparser2Reader():
                 # No elsewhere clauses were found so put the whole body into
                 # the single if block.
                 self.process_nodes(sched, node.content[1:-1])
-
         # Convert all uses of array syntax to indexed accesses
         self._array_syntax_to_indexed(ifblock, loop_vars)
         # Return the top-level loop generated by this handler


Test patch:
diff --git a/src/psyclone/tests/psyir/frontend/fparser2_test.py b/src/psyclone/tests/psyir/frontend/fparser2_test.py
index c04e6914d..932bda58c 100644
--- a/src/psyclone/tests/psyir/frontend/fparser2_test.py
+++ b/src/psyclone/tests/psyir/frontend/fparser2_test.py
@@ -58,8 +58,8 @@ from psyclone.psyir.frontend.fparser2 import (
 from psyclone.psyir.nodes import (
     Schedule, CodeBlock, Assignment, Return, UnaryOperation, BinaryOperation,
     IfBlock, Reference, ArrayReference, Literal, Range, KernelSchedule,
-    RegionDirective, Routine, StandaloneDirective, StructureReference,
-    ArrayOfStructuresReference, Call, IntrinsicCall)
+    RegionDirective, Routine, StandaloneDirective,
+    Call, IntrinsicCall)
 from psyclone.psyir.symbols import (
     DataSymbol, ContainerSymbol, ArgumentInterface, ArrayType,
     SymbolError, ScalarType, INTEGER_TYPE, REAL_TYPE, RoutineSymbol,
@@ -340,100 +340,6 @@ def test_get_arg_names(parser):
 # Class Fparser2Reader
 
 
-def test_array_notation_rank():
-    '''Test the static method _array_notation_rank in the fparser2reader
-    class.
-
-    '''
-    int_one = Literal("1", INTEGER_TYPE)
-    # Wrong type of argument
-    with pytest.raises(NotImplementedError) as err:
-        Fparser2Reader._array_notation_rank(int_one)
-    assert ("Expected either an ArrayReference, ArrayMember or a "
-            "StructureReference but got 'Literal'" in str(err.value))
-
-    # Structure reference containing no array access
-    symbol = DataSymbol("field", UnresolvedType())
-    with pytest.raises(InternalError) as err:
-        Fparser2Reader._array_notation_rank(
-            StructureReference.create(symbol, ["first", "second"]))
-    assert "No array access found in node 'field'" in str(err.value)
-
-    # Structure reference with ranges in more than one part reference.
-    lbound = IntrinsicCall.create(
-        IntrinsicCall.Intrinsic.LBOUND,
-        [StructureReference.create(symbol, ["first"]),
-         ("dim", int_one.copy())])
-    ubound = IntrinsicCall.create(
-        IntrinsicCall.Intrinsic.UBOUND,
-        [StructureReference.create(symbol, ["first"]),
-         ("dim", int_one.copy())])
-    my_range = Range.create(lbound, ubound)
-    with pytest.raises(InternalError) as err:
-        Fparser2Reader._array_notation_rank(
-            StructureReference.create(symbol, [("first", [my_range]),
-                                               ("second", [my_range.copy()])]))
-    assert ("Found a structure reference containing two or more part "
-            "references that have ranges: 'field%first(:)%second("
-            "LBOUND(field%first, dim=1):UBOUND(field%first, dim=1))'. This is "
-            "not valid within a WHERE in Fortran." in str(err.value))
-    # Repeat but this time for an ArrayOfStructuresReference.
-    with pytest.raises(InternalError) as err:
-        Fparser2Reader._array_notation_rank(
-            ArrayOfStructuresReference.create(symbol, [my_range.copy()],
-                                              ["first",
-                                               ("second", [my_range.copy()])]))
-    assert ("Found a structure reference containing two or more part "
-            "references that have ranges: 'field(LBOUND(field%first, dim=1):"
-            "UBOUND(field%first, dim=1))%first%second("
-            "LBOUND(field%first, dim=1):UBOUND(field%first, dim=1))'. This is "
-            "not valid within a WHERE in Fortran." in str(err.value))
-
-    # An array with no dimensions raises an exception
-    array_type = ArrayType(REAL_TYPE, [10])
-    symbol = DataSymbol("a", array_type)
-    array = ArrayReference(symbol)
-    with pytest.raises(InternalError) as excinfo:
-        Fparser2Reader._array_notation_rank(array)
-    assert ("ArrayReference malformed or incomplete: must have one or more "
-            "children representing array-index expressions but array 'a' has "
-            "none" in str(excinfo.value))
-
-    # If array syntax notation is found, it must be for all elements
-    # in that dimension
-    array_type = ArrayType(REAL_TYPE, [10, 10, 10])
-    symbol = DataSymbol("a", array_type)
-    lbound_op1 = IntrinsicCall.create(
-        IntrinsicCall.Intrinsic.LBOUND,
-        [Reference(symbol), ("dim", Literal("1", INTEGER_TYPE))])
-    ubound_op1 = IntrinsicCall.create(
-        IntrinsicCall.Intrinsic.UBOUND,
-        [Reference(symbol), ("dim", Literal("1", INTEGER_TYPE))])
-    lbound_op3 = IntrinsicCall.create(
-        IntrinsicCall.Intrinsic.LBOUND,
-        [Reference(symbol), ("dim", Literal("3", INTEGER_TYPE))])
-    ubound_op3 = IntrinsicCall.create(
-        IntrinsicCall.Intrinsic.UBOUND,
-        [Reference(symbol), ("dim", Literal("3", INTEGER_TYPE))])
-
-    range1 = Range.create(lbound_op1, ubound_op1)
-    range2 = Range.create(lbound_op3, ubound_op3)
-    one = Literal("1", INTEGER_TYPE)
-    array = ArrayReference.create(symbol, [range1, one.copy(), range2])
-    result = Fparser2Reader._array_notation_rank(array)
-    # Two array dimensions use array notation.
-    assert result == 2
-
-    # Make one of the array notation dimensions differ from what is required.
-    range2 = Range.create(lbound_op3.copy(), one.copy())
-    array = ArrayReference.create(symbol, [range1.copy(), one.copy(),
-                                           range2.copy()])
-    with pytest.raises(NotImplementedError) as excinfo:
-        Fparser2Reader._array_notation_rank(array)
-    assert ("Only array notation of the form my_array(:, :, ...) is "
-            "supported." in str(excinfo.value))
-
-
 def test_get_routine_schedules_wrong_module(parser):
     '''Test that get_routine_schedules() raises the expected errors if there
     are no or too many modules in the supplied parse tree.'''
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 a6ff7ca51..3b51510dc 100644
--- a/src/psyclone/tests/psyir/frontend/fparser2_where_handler_test.py
+++ b/src/psyclone/tests/psyir/frontend/fparser2_where_handler_test.py
@@ -45,12 +45,11 @@ from fparser.two import Fortran2003
 from psyclone.errors import InternalError
 from psyclone.psyir.frontend.fparser2 import Fparser2Reader
 from psyclone.psyir.nodes import (
-    ArrayMember, ArrayReference, Assignment, BinaryOperation, Call, CodeBlock,
-    Container, IfBlock, IntrinsicCall, Literal, Loop, Range, Routine, Schedule,
-    UnaryOperation)
+    ArrayMember, ArrayReference, Assignment, BinaryOperation,
+    Call, CodeBlock, Container, IfBlock, IntrinsicCall, Literal, Loop, Range,
+    Reference, Routine, Schedule, UnaryOperation)
 from psyclone.psyir.symbols import (
-    DataSymbol, ArrayType, ScalarType, REAL_TYPE, INTEGER_TYPE,
-    UnresolvedInterface)
+    DataSymbol, ScalarType, INTEGER_TYPE)
 
 
 def process_where(code, fparser_cls, symbols=None):
@@ -197,56 +196,6 @@ def test_labelled_where():
     assert isinstance(fake_parent.children[0], CodeBlock)
 
 
-@pytest.mark.usefixtures("parser")
-def test_missing_array_notation_lhs():
-    ''' Check that we get a code block if the WHERE does not use explicit
-    array syntax on the LHS of an assignment within the body.
-
-    '''
-    fake_parent, _ = process_where("WHERE (ptsu(:,:,:) /= 0._wp)\n"
-                                   "  z1_st = 1._wp / ptsu(:, :, :)\n"
-                                   "END WHERE\n", Fortran2003.Where_Construct,
-                                   ["ptsu", "z1_st"])
-    assert isinstance(fake_parent.children[0], CodeBlock)
-
-
-@pytest.mark.usefixtures("parser")
-def test_missing_array_notation_in_assign():
-    ''' Check that we get a code block if the WHERE contains an assignment
-    where no array notation appears. TODO #717 - extend this test so that
-    `z1_st` is of array type with rank 3. '''
-    fake_parent, _ = process_where("WHERE (ptsu(:,:,:) /= 0._wp)\n"
-                                   "  z1_st = 1._wp\n"
-                                   "END WHERE\n", Fortran2003.Where_Construct,
-                                   ["ptsu", "z1_st"])
-    assert isinstance(fake_parent.children[0], CodeBlock)
-
-
-@pytest.mark.usefixtures("parser")
-def test_where_array_notation_rank():
-    ''' Test that the _array_notation_rank() utility raises the expected
-    errors when passed an unsupported Array object.
-    '''
-    array_type = ArrayType(REAL_TYPE, [10])
-    symbol = DataSymbol("my_array", array_type)
-    my_array = ArrayReference(symbol)
-    processor = Fparser2Reader()
-    with pytest.raises(InternalError) as err:
-        processor._array_notation_rank(my_array)
-    assert ("ArrayReference malformed or incomplete: must have one or more "
-            "children representing array-index expressions but array "
-            "'my_array' has none." in str(err.value))
-    array_type = ArrayType(REAL_TYPE, [10])
-    my_array = ArrayReference.create(
-        DataSymbol("my_array", array_type),
-        [Range.create(Literal("2", INTEGER_TYPE),
-                      Literal("10", INTEGER_TYPE))])
-    with pytest.raises(NotImplementedError) as err:
-        processor._array_notation_rank(my_array)
-    assert ("Only array notation of the form my_array(:, :, ...) is "
-            "supported." in str(err.value))
-
-
 @pytest.mark.usefixtures("parser")
 def test_different_ranks_error():
     ''' Check that a WHERE construct containing array references of different
@@ -260,35 +209,6 @@ def test_different_ranks_error():
     assert isinstance(fake_parent.children[0], CodeBlock)
 
 
-@pytest.mark.usefixtures("parser")
-def test_array_notation_rank():
-    ''' Check that the _array_notation_rank() utility handles various examples
-    of array notation.
-
-    '''
-    fake_parent = Schedule()
-    fake_parent.symbol_table.new_symbol("z1_st")
-    fake_parent.symbol_table.new_symbol("ptsu")
-    fake_parent.symbol_table.new_symbol("n")
-    processor = Fparser2Reader()
-    reader = FortranStringReader("  z1_st(:, 2, :) = ptsu(:, :, 3)")
-    fparser2spec = Fortran2003.Assignment_Stmt(reader)
-    processor.process_nodes(fake_parent, [fparser2spec])
-    assert processor._array_notation_rank(fake_parent[0].lhs) == 2
-    reader = FortranStringReader("  z1_st(:, :, 2, :) = ptsu(:, :, :, 3)")
-    fparser2spec = Fortran2003.Assignment_Stmt(reader)
-    processor.process_nodes(fake_parent, [fparser2spec])
-    assert processor._array_notation_rank(fake_parent[1].lhs) == 3
-    # We don't support bounds on slices
-    reader = FortranStringReader("  z1_st(:, 1:n, 2, :) = ptsu(:, :, :, 3)")
-    fparser2spec = Fortran2003.Assignment_Stmt(reader)
-    processor.process_nodes(fake_parent, [fparser2spec])
-    with pytest.raises(NotImplementedError) as err:
-        processor._array_notation_rank(fake_parent[2].lhs)
-    assert ("Only array notation of the form my_array(:, :, ...) is "
-            "supported." in str(err.value))
-
-
 def test_where_symbol_clash(fortran_reader):
     ''' Check that we handle the case where the code we are processing
     already contains a symbol with the same name as one of the loop variables
@@ -318,8 +238,8 @@ def test_where_symbol_clash(fortran_reader):
 
 
 def test_where_within_loop(fortran_reader):
-    ''' Test for correct operation when we have a WHERE within an existing
-    loop and the referenced arrays are brought in from a module. '''
+    ''' Test for correct operation (Codeblock) when we have a WHERE within an
+    existing loop and the referenced arrays are brought in from a module. '''
     code = ("module my_mod\n"
             " use some_mod\n"
             "contains\n"
@@ -338,27 +258,8 @@ def test_where_within_loop(fortran_reader):
     mymod = psyir.children[0]
     assert isinstance(mymod, Container)
     sub = mymod.children[0]
-    assert "var" in sub.symbol_table
-    assert "var2" in sub.symbol_table
-    assert isinstance(sub.symbol_table.lookup("var").interface,
-                      UnresolvedInterface)
-    assert isinstance(sub.symbol_table.lookup("var2").interface,
-                      UnresolvedInterface)
     assert isinstance(sub, Routine)
-    assert isinstance(sub[0], Loop)
-    assert sub[0].variable.name == "jl"
-    where_loop = sub[0].loop_body[0]
-    assert isinstance(where_loop, Loop)
-    assert where_loop.variable.name == "widx1"
-    assert len(where_loop.loop_body.children) == 1
-    assert isinstance(where_loop.loop_body[0], IfBlock)
-    assign = where_loop.loop_body[0].if_body[0]
-    assert isinstance(assign, Assignment)
-    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)"
+    assert isinstance(sub[0].loop_body.children[0], CodeBlock)
 
 
 @pytest.mark.usefixtures("parser")
@@ -428,10 +329,17 @@ def test_where_mask_starting_value(fortran_reader, fortran_writer):
     code = '''\
     program my_sub
       use some_mod
+      type :: thing
+        real, dimension(100000, 100000, 1) :: z3
+      end type
+      integer, parameter :: jpl = 123435
+      integer, parameter :: jpr_ievp = 123
+      type(thing), dimension(10000) :: frcv
       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
+      integer, dimension(jpl) :: map
 
       WHERE( picefr(:,:) > 1.e-10 )
         zevap_ice(:,:,1) = snow(:,3,:) * frcv(jpr_ievp)%z3(:,:,1) / picefr(:,:)
@@ -447,8 +355,7 @@ def test_where_mask_starting_value(fortran_reader, fortran_writer):
     do widx1 = 1, 5 - (-5) + 1, 1
       if (picefr(-5 + widx1 - 1,-5 + widx2 - 1) > 1.e-10) then
         zevap_ice(widx1,widx2,1) = snow(-2 + widx1 - 1,3,-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) / \
+frcv(jpr_ievp)%z3(widx1,widx2,1) / \
 picefr(-5 + widx1 - 1,-5 + widx2 - 1)
       else
         zevap_ice(widx1,widx2,1) = snow(-2 + widx1 - 1,map(jpl),\
@@ -464,7 +371,14 @@ def test_where_mask_is_slice(fortran_reader, fortran_writer):
     '''
     code = '''\
     program my_sub
-      use some_mod
+      type :: thing
+        real, dimension(100000, 100000, 1) :: z3
+      end type
+      type(thing), dimension(10000) :: frcv
+      real, dimension(-5:5,-5:5) :: picefr
+      real, DIMENSION(11,11,jpl) :: zevap_ice
+      integer :: jstart, jstop, jpl, jpr_ievp
+
       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)
@@ -477,11 +391,8 @@ def test_where_mask_is_slice(fortran_reader, fortran_writer):
     # 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)
+    assert ("zevap_ice(widx1,widx2,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
@@ -495,8 +406,13 @@ def test_where_mask_is_slice_lower_limit(fortran_reader, fortran_writer):
     '''
     code = '''\
     subroutine my_sub(picefr)
-      use some_mod
       real, dimension(3:,2:) :: picefr
+      type :: thing
+        real, dimension(100000, 100000, 1) :: z3
+      end type
+      type(thing), dimension(10000) :: frcv
+      real, DIMENSION(11,11,jpl) :: zevap_ice
+      integer :: jstart, jstop, jpl, jpr_ievp
       WHERE( picefr(:4,jstart:) > 1.e-10 )
         zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) / picefr(:,:)
       ELSEWHERE ( picefr(4:5,jstart:) > 4.e-10)
@@ -513,8 +429,7 @@ def test_where_mask_is_slice_lower_limit(fortran_reader, fortran_writer):
     # specified in the original slice.
     assert ("if (picefr(LBOUND(picefr, dim=1) + 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)
+    assert ("zevap_ice(widx1,widx2,1) = 0.0" in out)
     assert "if (picefr(4 + widx1 - 1,jstart + widx2 - 1) > 4.e-10)" in out
 
 
@@ -563,6 +478,11 @@ def test_where_containing_sum_no_dim(fortran_reader, fortran_writer):
       REAL(wp), INTENT(in), DIMENSION(:,:) ::   picefr
       REAL(wp), DIMENSION(jpi,jpj,jpl) :: zevap_ice
       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple
+      type :: thing
+        real, dimension(100000, 100000, 1) :: z3
+      end type
+      integer :: jpr_ievp
+      type(thing), dimension(10000) :: frcv
       WHERE( picefr(:,:) > SUM(picefr) )
         zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) * &
                            SUM( a_i_last_couple ) / picefr(:,:)
@@ -604,19 +524,62 @@ def test_where_mask_containing_sum_with_dim(fortran_reader):
 
 
 @pytest.mark.usefixtures("parser")
-@pytest.mark.parametrize("rhs", ["depth", "maxval(depth(:))"])
-@pytest.mark.xfail(reason="#717 need to distinguish scalar and array "
-                   "assignments")
-def test_where_with_scalar_assignment(rhs):
+def test_where_with_scalar_assignment(fortran_reader, fortran_writer):
     ''' Test that a WHERE containing a scalar assignment is handled correctly.
-    Currently it is not as we do not distinguish between a scalar and an array
-    reference that is missing its colons. This will be fixed in #717.
+    '''
+    code = '''
+    subroutine sub()
+        integer, dimension(100,100,100) :: dry, z1_st, ptsu
+        integer :: var1, depth
+
+        where(dry(1, :, :))
+            var1 = depth
+            z1_st(:,2,:) = var1 / ptsu(:, :, 3)
+        end where
+    end subroutine sub
+    '''
+    fake_parent = fortran_reader.psyir_from_source(code)
+    # We should have a doubly-nested loop with an IfBlock inside
+    loops = fake_parent.walk(Loop)
+    assert len(loops) == 2
+    for loop in loops:
+        assert "was_where" in loop.annotations
+    assert isinstance(loops[1].loop_body[0], IfBlock)
+    out = fortran_writer(fake_parent)
+    correct = '''subroutine sub()
+  integer, dimension(100,100,100) :: dry
+  integer, dimension(100,100,100) :: z1_st
+  integer, dimension(100,100,100) :: ptsu
+  integer :: var1
+  integer :: depth
+  integer :: widx2
+  integer :: widx1
+
+  do widx2 = 1, SIZE(dry, dim=3), 1
+    do widx1 = 1, SIZE(dry, dim=2), 1
+      if (dry(1,widx1,widx2)) then
+        var1 = depth
+        z1_st(widx1,2,widx2) = var1 / ptsu(widx1,widx2,3)
+      end if
+    enddo
+  enddo
+
+end subroutine sub
+'''
+    assert out == correct
+
+
+@pytest.mark.usefixtures("parser")
+@pytest.mark.xfail(reason="#1960 Can't handle array-reduction intrinsics")
+def test_where_with_array_reduction_intrinsic():
+    ''' Test that a WHERE containing an array-reduction intrinsic is handled
+    correctly. Currently it is not supported. This will be fixed in #1960.
     '''
     fake_parent, _ = process_where(
-        f"WHERE (dry(1, :, :))\n"
-        f"  var1 = {rhs}\n"
-        f"  z1_st(:, 2, :) = var1 / ptsu(:, :, 3)\n"
-        f"END WHERE\n", Fortran2003.Where_Construct,
+        "WHERE (dry(1, :, :))\n"
+        "  var1 = maxval(depth(:))\n"
+        "  z1_st(:, 2, :) = var1 / ptsu(:, :, 3)\n"
+        "END WHERE\n", Fortran2003.Where_Construct,
         ["dry", "z1_st", "depth", "ptsu", "var1"])
     # We should have a doubly-nested loop with an IfBlock inside
     loops = fake_parent.walk(Loop)
@@ -800,19 +763,34 @@ def test_where_ordering(parser):
     [("where (my_type%var(:) > epsi20)\n"
       "my_type%array(:,jl) = 3.0\n", "my_type%var"),
      ("where (my_type%var(:) > epsi20)\n"
-      "my_type(jl)%array(:,jl) = 3.0\n", "my_type%var"),
+      "my_type2(jl)%array2(:,jl) = 3.0\n", "my_type%var"),
      ("where (my_type%block(jl)%var(:) > epsi20)\n"
       "my_type%block%array(:,jl) = 3.0\n", "my_type%block(jl)%var"),
      ("where (my_type%block(jl)%var(:) > epsi20)\n"
-      "my_type%block(jl)%array(:,jl) = 3.0\n", "my_type%block(jl)%var")])
-def test_where_derived_type(fortran_reader, fortran_writer, code, size_arg):
+      "my_type%block(jl)%array(:,jl) = 3.0\n", "my_type%block(jl)%var"),
+     ("where (my_type2(:)%var > epsi20)\n"
+      "my_type2(:)%var = 3.0\n", "my_type2")])
+def test_where_derived_type(fortran_reader, code, size_arg):
     ''' Test that we handle the case where array members of a derived type
     are accessed within a WHERE. '''
     code = (f"module my_mod\n"
             f" use some_mod\n"
             f"contains\n"
             f"subroutine my_sub()\n"
+            f"  type :: subtype\n"
+            f"    real, dimension(:) :: var\n"
+            f"    real, dimension(:,:) :: array\n"
+            f" end type\n"
+            f"  type :: sometype\n"
+            f"    real, dimension(:) :: var\n"
+            f"    real, dimension(:) :: array\n"
+            f"    real, dimension(:,:) :: array2\n"
+            f"    type(subtype), dimension(:) :: block\n"
+            f"  end type\n"
+            f"  type(sometype) :: my_type\n"
+            f"  type(sometype), dimension(:) :: my_type2\n"
             f"  integer :: jl\n"
+            f"  real :: epsi20\n"
             f"  do jl = 1, 10\n"
             f"{code}"
             f"    end where\n"
@@ -831,3 +809,198 @@ def test_where_derived_type(fortran_reader, fortran_writer, code, size_arg):
     array_members = loops[0].walk(ArrayMember)
     for member in array_members:
         assert "+ widx1 - 1" in member.indices[0].debug_string()
+
+
+@pytest.mark.parametrize(
+    "code",
+    [("where (my_type%var > epsi20)\n"
+      "my_type%array(:,jl) = 3.0\n"),
+     ("where (my_type%var > epsi20)\n"
+      "my_type(jl)%array(:,jl) = 3.0\n"),
+     ("where (my_type%block(jl)%var > epsi20)\n"
+      "my_type(jl%array(:,jl) = 3.0\n"),
+     ("where (my_type%block(jl)%var > epsi20)\n"
+      "my_type%block(jl)%var = 3.0\n")])
+@pytest.mark.xfail(reason="#1960 Can't handle WHERE constructs without "
+                          "explicit array notation inside derived types.")
+def test_where_noarray_syntax_derived_types(fortran_reader, code):
+    '''Xfailing test for when a derived type access in a where condition
+    doesn't use range syntax.'''
+    code = (f"module my_mod\n"
+            f" use some_mod\n"
+            f"contains\n"
+            f"subroutine my_sub()\n"
+            f"  integer :: jl\n"
+            f"  do jl = 1, 10\n"
+            f"{code}"
+            f"    end where\n"
+            f"  end do\n"
+            f"end subroutine my_sub\n"
+            f"end module my_mod\n")
+    psyir = fortran_reader.psyir_from_source(code)
+    cbs = psyir.walk(CodeBlock)
+    assert len(cbs) == 0
+
+
+@pytest.mark.parametrize(
+    "code",
+    [("where (var(:) > 1.0)\n"
+        "var(:) = 3.0\n"),
+     ("where (var(:) > var2(1:20))\n"
+      "var(:) = 3.0\n"),
+     ("where (var > 1.0)\n"
+      "var = 3.0\n"),
+     ("where (var > var2)\n"
+      "var = 3.0\n"),
+     ("where (var(:) > var2)\n"
+      "var(:) = 3.0\n"),
+     ("where (var > var2(1:20))\n"
+      "var = 3.0\n")])
+def test_where_scalar_var(fortran_reader, code):
+    '''Test where we have a scalar variable in a WHERE clause with no
+    array index clause.'''
+    code = (
+        f"program where_test\n"
+        f"implicit none\n"
+        f"integer, parameter :: nc=20\n"
+        f"integer :: ii\n"
+        f"real :: var(nc)\n"
+        f"real:: var2(nc)\n"
+        f"var = 1.5\n"
+        f"var2 = 1.6\n"
+        f"{code}"
+        f"end where\n"
+        f"end program where_test")
+    psyir = fortran_reader.psyir_from_source(code)
+    loops = psyir.walk(Loop)
+    assert len(loops) == 1
+    assert isinstance(loops[0].stop_expr, Reference)
+    assert loops[0].stop_expr.debug_string() == "nc"
+    assert isinstance(loops[0].loop_body[0], IfBlock)
+    # All Range nodes should have been replaced
+    assert not loops[0].walk(Range)
+    array_refs = loops[0].walk(ArrayReference)
+    for member in array_refs:
+        assert "widx1" in member.indices[0].debug_string()
+
+    # All References to var or var2 should be ArrayReferences
+    refs = loops[0].walk(Reference)
+    for ref in refs:
+        if "var" in ref.symbol.name:
+            assert isinstance(ref, ArrayReference)
+
+
+def test_import_in_where_clause(fortran_reader):
+    '''
+    Test that imported symbols inside a where clause produce a code block
+    as we don't know what is scalar vs an array.
+    '''
+    code = '''
+    program where_test
+    implicit none
+    use some_module, only: a, b, c, d
+
+    where(a(:) + b > 1)
+       b = c + d
+    end where
+    end program
+    '''
+    psyir = fortran_reader.psyir_from_source(code)
+    assert isinstance(psyir.children[0].children[0], CodeBlock)
+
+    code2 = '''
+    program where_test
+    implicit none
+    use some_module, only: c, d
+    integer, dimension(100) :: a, b
+
+    where(a(:) + b > 1)
+       b = c + d
+    end where
+    end program
+    '''
+    psyir = fortran_reader.psyir_from_source(code2)
+    assert isinstance(psyir.children[0].children[0], CodeBlock)
+
+
+def test_non_array_reduction_intrinsic(fortran_reader, fortran_writer):
+    '''
+    Test that a non-array reduction intrinsic (such as SIN) produces
+    the correct PSyIR.
+    '''
+    code = '''
+    program where_test
+    implicit none
+    integer, dimension(100) :: a, b, c
+
+    where(a < b)
+       c = sin(a)
+    end where
+    end program
+    '''
+    psyir = fortran_reader.psyir_from_source(code)
+    assert isinstance(psyir.children[0].children[0], Loop)
+    out = fortran_writer(psyir)
+    correct = '''program where_test
+  integer, dimension(100) :: a
+  integer, dimension(100) :: b
+  integer, dimension(100) :: c
+  integer :: widx1
+
+  do widx1 = 1, 100, 1
+    if (a(widx1) < b(widx1)) then
+      c(widx1) = SIN(a(widx1))
+    end if
+  enddo
+
+end program where_test'''
+    assert correct in out
+
+
+def test_non_elemental_intrinsic(fortran_reader):
+    '''
+    Test that a non-elemental reduction intrinsic (such as DOT_PRODUCT)
+    produces the correct PSyIR.
+    '''
+    code = '''
+    program where_test
+    implicit none
+    integer, dimension(100) :: a, b, c
+
+    where(a < b)
+       c = DOT_PRODUCT(a, b)
+    end where
+    end program
+    '''
+    psyir = fortran_reader.psyir_from_source(code)
+    # This should be a Loop I suppose, this is valid Fortran
+    # as long as a and b aren't expanded
+    assert isinstance(psyir.children[0].children[0], Loop)
+    intrinsic = psyir.walk(IntrinsicCall)[0]
+    assert intrinsic.intrinsic.name == "DOT_PRODUCT"
+    assert intrinsic.children[1].name == "a"
+    assert not isinstance(intrinsic.children[1], ArrayReference)
+    assert intrinsic.children[2].name == "b"
+    assert not isinstance(intrinsic.children[2], ArrayReference)
+
+
+def test_non_array_ref_intrinsic_transformation_error(fortran_reader):
+    '''
+    Test for coverage for the try/except in the intrinsic_ancestor section
+    of _array_syntax_to_indexed sub function.'''
+
+    code = '''
+    program where_test
+    implicit none
+    integer, dimension(100) :: a, b, c
+    real :: d
+
+    where(a < b)
+       c = sin(d)
+    end where
+    end program
+    '''
+    psyir = fortran_reader.psyir_from_source(code)
+    references = psyir.walk(Reference)
+    # The d should not have been transformed into an array.
+    assert not isinstance(references[7], ArrayReference)