diff --git a/flang/include/flang/Evaluate/constant.h b/flang/include/flang/Evaluate/constant.h
--- a/flang/include/flang/Evaluate/constant.h
+++ b/flang/include/flang/Evaluate/constant.h
@@ -140,7 +140,8 @@
     }
   }
 
-  // Apply subscripts.
+  // Apply subscripts.  An empty subscript list is allowed for
+  // a scalar constant.
   Element At(const ConstantSubscripts &) const;
 
   Constant Reshape(ConstantSubscripts &&) const;
@@ -177,7 +178,7 @@
     }
   }
 
-  // Apply subscripts
+  // Apply subscripts, if any.
   Scalar<Result> At(const ConstantSubscripts &) const;
 
   Constant Reshape(ConstantSubscripts &&) const;
diff --git a/flang/lib/Evaluate/fold-implementation.h b/flang/lib/Evaluate/fold-implementation.h
--- a/flang/lib/Evaluate/fold-implementation.h
+++ b/flang/lib/Evaluate/fold-implementation.h
@@ -443,9 +443,8 @@
     // Compute the shape of the result based on shapes of arguments
     ConstantSubscripts shape;
     int rank{0};
-    const ConstantSubscripts *shapes[sizeof...(TA)]{
-        &std::get<I>(*args)->shape()...};
-    const int ranks[sizeof...(TA)]{std::get<I>(*args)->Rank()...};
+    const ConstantSubscripts *shapes[]{&std::get<I>(*args)->shape()...};
+    const int ranks[]{std::get<I>(*args)->Rank()...};
     for (unsigned int i{0}; i < sizeof...(TA); ++i) {
       if (ranks[i] > 0) {
         if (rank == 0) {
@@ -470,20 +469,19 @@
     std::vector<Scalar<TR>> results;
     if (TotalElementCount(shape) > 0) {
       ConstantBounds bounds{shape};
-      ConstantSubscripts index(rank, 1);
+      ConstantSubscripts resultIndex(rank, 1);
+      ConstantSubscripts argIndex[]{std::get<I>(*args)->lbounds()...};
       do {
         if constexpr (std::is_same_v<WrapperType<TR, TA...>,
                           ScalarFuncWithContext<TR, TA...>>) {
-          results.emplace_back(func(context,
-              (ranks[I] ? std::get<I>(*args)->At(index)
-                        : std::get<I>(*args)->GetScalarValue().value())...));
+          results.emplace_back(
+              func(context, std::get<I>(*args)->At(argIndex[I])...));
         } else if constexpr (std::is_same_v<WrapperType<TR, TA...>,
                                  ScalarFunc<TR, TA...>>) {
-          results.emplace_back(func(
-              (ranks[I] ? std::get<I>(*args)->At(index)
-                        : std::get<I>(*args)->GetScalarValue().value())...));
+          results.emplace_back(func(std::get<I>(*args)->At(argIndex[I])...));
         }
-      } while (bounds.IncrementSubscripts(index));
+        (std::get<I>(*args)->IncrementSubscripts(argIndex[I]), ...);
+      } while (bounds.IncrementSubscripts(resultIndex));
     }
     // Build and return constant result
     if constexpr (TR::category == TypeCategory::Character) {
@@ -732,17 +730,11 @@
     Expr<T> folded{Fold(context_, common::Clone(expr.value()))};
     if (const auto *c{UnwrapConstantValue<T>(folded)}) {
       // Copy elements in Fortran array element order
-      ConstantSubscripts shape{c->shape()};
-      int rank{c->Rank()};
-      ConstantSubscripts index(GetRank(shape), 1);
-      for (std::size_t n{c->size()}; n-- > 0;) {
-        elements_.emplace_back(c->At(index));
-        for (int d{0}; d < rank; ++d) {
-          if (++index[d] <= shape[d]) {
-            break;
-          }
-          index[d] = 1;
-        }
+      if (c->size() > 0) {
+        ConstantSubscripts index{c->lbounds()};
+        do {
+          elements_.emplace_back(c->At(index));
+        } while (c->IncrementSubscripts(index));
       }
       return true;
     } else {
diff --git a/flang/test/Evaluate/folding16.f90 b/flang/test/Evaluate/folding16.f90
new file mode 100644
--- /dev/null
+++ b/flang/test/Evaluate/folding16.f90
@@ -0,0 +1,8 @@
+! RUN: %S/test_folding.sh %s %t %f18
+! Ensure that lower bounds are accounted for in intrinsic folding;
+! this is a regression test for a bug in which they were not
+real, parameter :: a(-1:-1) = 1.
+real, parameter :: b(-1:-1) = log(a)
+logical, parameter :: test = lbound(a,1)==-1 .and. lbound(b,1)==-1 .and. &
+                             lbound(log(a),1)==1 .and. all(b==0)
+end