diff --git a/changelog b/changelog
index 3f43a5e0c0..793ae1557a 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,6 @@
+   6) PR #3369 towards #1658. Switches NEMO scripts to use all of the
+   supported reduction intrinsic expansion transformations.
+
    5) PR #3377 for #2612. Add support for Fortran character length.
 
    4) PR #3388 for #3334. Remove the functionality to write (and rename)
diff --git a/examples/nemo/scripts/utils.py b/examples/nemo/scripts/utils.py
index 3cab05f169..fec2f529f5 100755
--- a/examples/nemo/scripts/utils.py
+++ b/examples/nemo/scripts/utils.py
@@ -45,9 +45,10 @@
 from psyclone.psyir.symbols import DataSymbol
 from psyclone.psyir.transformations import (
     ArrayAssignment2LoopsTrans, HoistLoopBoundExprTrans, HoistLocalArraysTrans,
-    HoistTrans, InlineTrans, Maxval2LoopTrans, ProfileTrans,
-    OMPMinimiseSyncTrans, Reference2ArrayRangeTrans,
-    ScalarisationTrans, IncreaseRankLoopArraysTrans, MaximalRegionTrans)
+    HoistTrans, InlineTrans, Maxval2LoopTrans, Sum2LoopTrans, Minval2LoopTrans,
+    Product2LoopTrans, ProfileTrans, OMPMinimiseSyncTrans,
+    Reference2ArrayRangeTrans, ScalarisationTrans, IncreaseRankLoopArraysTrans,
+    MaximalRegionTrans)
 from psyclone.transformations import TransformationError
 
 # USE statements to chase to gather additional symbol information.
@@ -219,11 +220,17 @@ def normalise_loops(
 
     if loopify_array_intrinsics:
         for intr in schedule.walk(IntrinsicCall):
-            if intr.intrinsic.name == "MAXVAL":
-                try:
+            try:
+                if intr.intrinsic.name == "MAXVAL":
                     Maxval2LoopTrans().apply(intr, verbose=True)
-                except TransformationError as err:
-                    print(err.value)
+                elif intr.intrinsic.name == "SUM":
+                    Sum2LoopTrans().apply(intr, verbose=True)
+                elif intr.intrinsic.name == "MINVAL":
+                    Minval2LoopTrans().apply(intr, verbose=True)
+                elif intr.intrinsic.name == "PRODUCT":
+                    Product2LoopTrans().apply(intr, verbose=True)
+            except TransformationError as err:
+                print(err.value)
 
     if convert_range_loops:
         # Convert all array implicit loops to explicit loops
diff --git a/src/psyclone/psyir/nodes/intrinsic_call.py b/src/psyclone/psyir/nodes/intrinsic_call.py
index 9b132ca149..9b42a183b2 100644
--- a/src/psyclone/psyir/nodes/intrinsic_call.py
+++ b/src/psyclone/psyir/nodes/intrinsic_call.py
@@ -239,28 +239,29 @@ def _type_of_named_arg_with_optional_kind_and_dim(
     return _type_of_arg_with_rank_minus_one(arg, dtype)
 
 
-def _type_with_specified_precision_and_optional_dim(
+def _type_of_named_arg_accounting_for_dim_arg(
         node: IntrinsicCall, argument_name: str,
-        intrinsic: ScalarType.Intrinsic = ScalarType.Intrinsic.BOOLEAN
-        ) -> DataType:
+) -> DataType:
     """Helper function for the common IntrinsicCall case where the
-    return type is a Scalar with the precision of a named argument,
-    unless an optional argument named 'dim' exists, in which case an array
-    with rank one less than the input node is given instead.
+    return type is the same as the given named argument. If intrinsiccall
+    has no 'dim' argument, it returns its elemental type, but if 'dim'
+    exists, it will be the given named argument's rank minus one.
 
-    :param node: The IntrinsicCall whose return type to compute.
-    :param argument_name: The name of the argument whose precision to be used.
-    :param intrinsic: The type of the intrinsic of the resulting datatype.
-                      Default is ScalarType.Intrinsic.BOOLEAN
+    :param node: the IntrinsicCall whose return type to compute.
+    :param argument_name: the name of the argument whose type to use.
 
     :returns: the computed datatype for the IntrinsicCall.
     """
-    dtype = ScalarType(
-        intrinsic, node.argument_by_name(argument_name).datatype.precision
-    )
-    # If dim is not present, or the rank of the
-    # array argument is 1 then this returns a scalar.
     arg = node.argument_by_name(argument_name)
+    arg_dt = arg.datatype
+    if (
+        not isinstance(arg_dt, ArrayType) or
+        not isinstance(arg_dt.elemental_type, ScalarType) or
+        not isinstance(arg_dt.elemental_type.intrinsic, ScalarType.Intrinsic)
+    ):
+        return UnresolvedType()
+    dtype = arg_dt.elemental_type
+    # If dim is not present, return the same datatype
     if "dim" not in node.argument_names:
         return dtype
 
@@ -877,7 +878,7 @@ class Intrinsic(IAttr, Enum):
             optional_args={"dim": DataNode},
             return_type=(
                 lambda node:
-                _type_with_specified_precision_and_optional_dim(
+                _type_of_named_arg_accounting_for_dim_arg(
                     node, "mask"
                 )
             ),
@@ -945,7 +946,7 @@ class Intrinsic(IAttr, Enum):
             optional_args={"dim": DataNode},
             return_type=(
                 lambda node:
-                _type_with_specified_precision_and_optional_dim(
+                _type_of_named_arg_accounting_for_dim_arg(
                     node, "mask"
                 )
             ),
@@ -3885,9 +3886,8 @@ class Intrinsic(IAttr, Enum):
             optional_args={"mask": DataNode},
             return_type=(
                 lambda node:
-                _type_with_specified_precision_and_optional_dim(
+                _type_of_named_arg_accounting_for_dim_arg(
                     node, "array",
-                    node.argument_by_name("array").datatype.intrinsic
                 )
             ),
             reference_accesses=lambda node: (
@@ -4542,9 +4542,8 @@ class Intrinsic(IAttr, Enum):
             optional_args={"mask": DataNode},
             return_type=(
                 lambda node:
-                _type_with_specified_precision_and_optional_dim(
+                _type_of_named_arg_accounting_for_dim_arg(
                     node, "array",
-                    node.argument_by_name("array").datatype.intrinsic
                 )
             ),
             reference_accesses=lambda node: (
diff --git a/src/psyclone/psyir/symbols/symbol.py b/src/psyclone/psyir/symbols/symbol.py
index be8f863f8b..1580b446a7 100644
--- a/src/psyclone/psyir/symbols/symbol.py
+++ b/src/psyclone/psyir/symbols/symbol.py
@@ -499,7 +499,7 @@ def is_array_access(self, index_variable=None, access_info=None):
             access information is given.
 
         '''
-        # TODO #1270: this function might either be better off elsewhere,
+        # TODO #3098: this function might either be better off elsewhere,
         # or even do not implement one function that uses both access
         # information and symbol table - if required, the user can
         # query both in two simple statements anyway.
diff --git a/src/psyclone/psyir/tools/dependency_tools.py b/src/psyclone/psyir/tools/dependency_tools.py
index 7cf69bcbfb..5bc8a3d9b4 100644
--- a/src/psyclone/psyir/tools/dependency_tools.py
+++ b/src/psyclone/psyir/tools/dependency_tools.py
@@ -898,9 +898,11 @@ def can_loop_be_parallelised(self, loop,
             except AttributeError:
                 # If its a node without a symbol, look it up
                 var_name = signature.var_name
-                symbol = symbol_table.lookup(var_name)
+                symbol = symbol_table.lookup(var_name, otherwise=None)
+                if symbol is None:
+                    return False
 
-            # TODO #1270 - the is_array_access function might be moved
+            # TODO #3098 - the is_array_access function might be moved
             is_array = symbol.is_array_access(access_info=var_info)
             if is_array:
                 # Handle arrays
@@ -983,7 +985,7 @@ def can_loops_be_fused(self, loop1, loop2):
                 continue
 
             symbol = symbol_table.lookup(signature.var_name)
-            # TODO #1270 - the is_array_access function might be moved
+            # TODO #3098 - the is_array_access function might be moved
             is_array = symbol.is_array_access(access_info=var_info1)
             if not is_array:
                 result = self._fuse_validate_written_scalar(var_info1,
diff --git a/src/psyclone/tests/psyir/nodes/intrinsic_call_test.py b/src/psyclone/tests/psyir/nodes/intrinsic_call_test.py
index 8e5a4f232f..9b760c2ed6 100644
--- a/src/psyclone/tests/psyir/nodes/intrinsic_call_test.py
+++ b/src/psyclone/tests/psyir/nodes/intrinsic_call_test.py
@@ -61,7 +61,7 @@
     _type_of_arg_with_rank_minus_one,
     _type_of_named_argument,
     _type_of_named_arg_with_optional_kind_and_dim,
-    _type_with_specified_precision_and_optional_dim,
+    _type_of_named_arg_accounting_for_dim_arg,
     _type_of_scalar_with_optional_kind,
     _type_of_intrinsic_with_argname_kind_and_optional_dim,
     _type_of_intrinsic_with_precision_of_named_arg,
@@ -1064,27 +1064,29 @@ def test_type_of_named_arg_with_optional_kind_and_dim(
     assert dtype.precision.value == "8"
 
 
-def test_type_with_specified_precision_and_optional_dim(fortran_reader):
-    """Test the _type_with_specified_precision_and_optional_dim
+def test_type_of_named_arg_accounting_for_dim_arg(fortran_reader):
+    """Test the _type_of_named_arg_accounting_for_dim_arg
     helper function."""
     code = """subroutine test
+    use other
     integer, dimension(100, 100) :: x
     integer :: y
     y = PRODUCT(x)
     y = PRODUCT(x, dim=2)
+    y = PRODUCT(z, dim=2)
     end subroutine test"""
     psyir = fortran_reader.psyir_from_source(code)
     intrinsics = psyir.walk(IntrinsicCall)
 
-    dtype = _type_with_specified_precision_and_optional_dim(
-        intrinsics[0], "array", ScalarType.Intrinsic.INTEGER,
+    dtype = _type_of_named_arg_accounting_for_dim_arg(
+        intrinsics[0], "array"
     )
     assert isinstance(dtype, ScalarType)
     assert dtype.intrinsic == ScalarType.Intrinsic.INTEGER
     assert dtype.precision == ScalarType.Precision.UNDEFINED
 
-    dtype = _type_with_specified_precision_and_optional_dim(
-        intrinsics[1], "array", ScalarType.Intrinsic.INTEGER,
+    dtype = _type_of_named_arg_accounting_for_dim_arg(
+        intrinsics[1], "array"
     )
     assert isinstance(dtype, ArrayType)
     assert len(dtype.shape) == 1
@@ -1092,6 +1094,11 @@ def test_type_with_specified_precision_and_optional_dim(fortran_reader):
     assert dtype.intrinsic == ScalarType.Intrinsic.INTEGER
     assert dtype.precision == ScalarType.Precision.UNDEFINED
 
+    dtype = _type_of_named_arg_accounting_for_dim_arg(
+        intrinsics[2], "array"
+    )
+    assert isinstance(dtype, UnresolvedType)
+
 
 def test_type_of_intrinsic_with_precision_of_named_arg(fortran_reader):
     """Test the _type_of_intrinsic_with_precision_of_named_arg helper