diff --git a/flang-rt/lib/runtime/numeric.cpp b/flang-rt/lib/runtime/numeric.cpp
index 37638765dc65..78f148dbc5d8 100644
--- a/flang-rt/lib/runtime/numeric.cpp
+++ b/flang-rt/lib/runtime/numeric.cpp
@@ -229,6 +229,24 @@ RT_API_ATTRS BTy FPowI(BTy base, ETy exp) {
   return result;
 }
 
+// Exponentiation operator for (Unsigned ** Unsigned) cases
+template <typename Ty> RT_API_ATTRS Ty UPow(Ty base, Ty exp) {
+  if (exp == Ty{0})
+    return Ty{1};
+  Ty result{1};
+  while (true) {
+    if (exp & Ty{1}) {
+      result *= base;
+    }
+    exp >>= 1;
+    if (exp == Ty{0}) {
+      break;
+    }
+    base *= base;
+  }
+  return result;
+}
+
 extern "C" {
 RT_EXT_API_GROUP_BEGIN
 
@@ -933,6 +951,27 @@ CppTypeFor<TypeCategory::Real, 16> RTDEF(FPow16k)(
 }
 #endif
 
+CppTypeFor<TypeCategory::Unsigned, 1> RTDEF(UPow1)(
+    CppTypeFor<TypeCategory::Unsigned, 1> b,
+    CppTypeFor<TypeCategory::Unsigned, 1> e) {
+  return UPow(b, e);
+}
+CppTypeFor<TypeCategory::Unsigned, 2> RTDEF(UPow2)(
+    CppTypeFor<TypeCategory::Unsigned, 2> b,
+    CppTypeFor<TypeCategory::Unsigned, 2> e) {
+  return UPow(b, e);
+}
+CppTypeFor<TypeCategory::Unsigned, 4> RTDEF(UPow4)(
+    CppTypeFor<TypeCategory::Unsigned, 4> b,
+    CppTypeFor<TypeCategory::Unsigned, 4> e) {
+  return UPow(b, e);
+}
+CppTypeFor<TypeCategory::Unsigned, 8> RTDEF(UPow8)(
+    CppTypeFor<TypeCategory::Unsigned, 8> b,
+    CppTypeFor<TypeCategory::Unsigned, 8> e) {
+  return UPow(b, e);
+}
+
 RT_EXT_API_GROUP_END
 } // extern "C"
 } // namespace Fortran::runtime
diff --git a/flang/include/flang/Evaluate/expression.h b/flang/include/flang/Evaluate/expression.h
index 1203fca8640a..f7a1f9b95518 100644
--- a/flang/include/flang/Evaluate/expression.h
+++ b/flang/include/flang/Evaluate/expression.h
@@ -566,9 +566,9 @@ private:
   using Conversions = std::tuple<Convert<Result, TypeCategory::Integer>,
       Convert<Result, TypeCategory::Real>,
       Convert<Result, TypeCategory::Unsigned>>;
-  using Operations =
-      std::tuple<Parentheses<Result>, Negate<Result>, Add<Result>,
-          Subtract<Result>, Multiply<Result>, Divide<Result>, Extremum<Result>>;
+  using Operations = std::tuple<Parentheses<Result>, Negate<Result>,
+      Add<Result>, Subtract<Result>, Multiply<Result>, Divide<Result>,
+      Power<Result>, Extremum<Result>>;
   using Others = std::tuple<Constant<Result>, ArrayConstructor<Result>,
       Designator<Result>, FunctionRef<Result>>;
 
diff --git a/flang/include/flang/Evaluate/tools.h b/flang/include/flang/Evaluate/tools.h
index 7fcd04778301..225e1a77f041 100644
--- a/flang/include/flang/Evaluate/tools.h
+++ b/flang/include/flang/Evaluate/tools.h
@@ -750,11 +750,11 @@ Expr<SomeKind<CAT>> PromoteAndCombine(
 // one of the operands to the type of the other.  Handles special cases with
 // typeless literal operands and with REAL/COMPLEX exponentiation to INTEGER
 // powers.
-template <template <typename> class OPR, bool CAN_BE_UNSIGNED = true>
+template <template <typename> class OPR>
 std::optional<Expr<SomeType>> NumericOperation(parser::ContextualMessages &,
     Expr<SomeType> &&, Expr<SomeType> &&, int defaultRealKind);
 
-extern template std::optional<Expr<SomeType>> NumericOperation<Power, false>(
+extern template std::optional<Expr<SomeType>> NumericOperation<Power>(
     parser::ContextualMessages &, Expr<SomeType> &&, Expr<SomeType> &&,
     int defaultRealKind);
 extern template std::optional<Expr<SomeType>> NumericOperation<Multiply>(
diff --git a/flang/include/flang/Runtime/numeric.h b/flang/include/flang/Runtime/numeric.h
index 794c8f476882..17ed31ab8790 100644
--- a/flang/include/flang/Runtime/numeric.h
+++ b/flang/include/flang/Runtime/numeric.h
@@ -453,6 +453,19 @@ CppTypeFor<TypeCategory::Real, 16> RTDECL(FPow16k)(
     CppTypeFor<TypeCategory::Integer, 8> e);
 #endif
 
+CppTypeFor<TypeCategory::Unsigned, 1> RTDEF(UPow1)(
+    CppTypeFor<TypeCategory::Unsigned, 1> b,
+    CppTypeFor<TypeCategory::Unsigned, 1> e);
+CppTypeFor<TypeCategory::Unsigned, 2> RTDEF(UPow2)(
+    CppTypeFor<TypeCategory::Unsigned, 2> b,
+    CppTypeFor<TypeCategory::Unsigned, 2> e);
+CppTypeFor<TypeCategory::Unsigned, 4> RTDEF(UPow4)(
+    CppTypeFor<TypeCategory::Unsigned, 4> b,
+    CppTypeFor<TypeCategory::Unsigned, 4> e);
+CppTypeFor<TypeCategory::Unsigned, 8> RTDEF(UPow8)(
+    CppTypeFor<TypeCategory::Unsigned, 8> b,
+    CppTypeFor<TypeCategory::Unsigned, 8> e);
+
 } // extern "C"
 } // namespace Fortran::runtime
 #endif // FORTRAN_RUNTIME_NUMERIC_H_
diff --git a/flang/lib/Evaluate/tools.cpp b/flang/lib/Evaluate/tools.cpp
index e6e215ecb489..1f3cbbf6a0c3 100644
--- a/flang/lib/Evaluate/tools.cpp
+++ b/flang/lib/Evaluate/tools.cpp
@@ -495,7 +495,7 @@ Expr<SomeComplex> PromoteMixedComplexReal(
 // N.B. When a "typeless" BOZ literal constant appears as one (not both!) of
 // the operands to a dyadic operation where one is permitted, it assumes the
 // type and kind of the other operand.
-template <template <typename> class OPR, bool CAN_BE_UNSIGNED>
+template <template <typename> class OPR>
 std::optional<Expr<SomeType>> NumericOperation(
     parser::ContextualMessages &messages, Expr<SomeType> &&x,
     Expr<SomeType> &&y, int defaultRealKind) {
@@ -510,13 +510,8 @@ std::optional<Expr<SomeType>> NumericOperation(
                 std::move(rx), std::move(ry)));
           },
           [&](Expr<SomeUnsigned> &&ix, Expr<SomeUnsigned> &&iy) {
-            if constexpr (CAN_BE_UNSIGNED) {
-              return Package(PromoteAndCombine<OPR, TypeCategory::Unsigned>(
-                  std::move(ix), std::move(iy)));
-            } else {
-              messages.Say("Operands must not be UNSIGNED"_err_en_US);
-              return NoExpr();
-            }
+            return Package(PromoteAndCombine<OPR, TypeCategory::Unsigned>(
+                std::move(ix), std::move(iy)));
           },
           // Mixed REAL/INTEGER operations
           [](Expr<SomeReal> &&rx, Expr<SomeInteger> &&iy) {
@@ -575,34 +570,31 @@ std::optional<Expr<SomeType>> NumericOperation(
           },
           // Operations with one typeless operand
           [&](BOZLiteralConstant &&bx, Expr<SomeInteger> &&iy) {
-            return NumericOperation<OPR, CAN_BE_UNSIGNED>(messages,
+            return NumericOperation<OPR>(messages,
                 AsGenericExpr(ConvertTo(iy, std::move(bx))), std::move(y),
                 defaultRealKind);
           },
           [&](BOZLiteralConstant &&bx, Expr<SomeUnsigned> &&iy) {
-            return NumericOperation<OPR, CAN_BE_UNSIGNED>(messages,
+            return NumericOperation<OPR>(messages,
                 AsGenericExpr(ConvertTo(iy, std::move(bx))), std::move(y),
                 defaultRealKind);
           },
           [&](BOZLiteralConstant &&bx, Expr<SomeReal> &&ry) {
-            return NumericOperation<OPR, CAN_BE_UNSIGNED>(messages,
+            return NumericOperation<OPR>(messages,
                 AsGenericExpr(ConvertTo(ry, std::move(bx))), std::move(y),
                 defaultRealKind);
           },
           [&](Expr<SomeInteger> &&ix, BOZLiteralConstant &&by) {
-            return NumericOperation<OPR, CAN_BE_UNSIGNED>(messages,
-                std::move(x), AsGenericExpr(ConvertTo(ix, std::move(by))),
-                defaultRealKind);
+            return NumericOperation<OPR>(messages, std::move(x),
+                AsGenericExpr(ConvertTo(ix, std::move(by))), defaultRealKind);
           },
           [&](Expr<SomeUnsigned> &&ix, BOZLiteralConstant &&by) {
-            return NumericOperation<OPR, CAN_BE_UNSIGNED>(messages,
-                std::move(x), AsGenericExpr(ConvertTo(ix, std::move(by))),
-                defaultRealKind);
+            return NumericOperation<OPR>(messages, std::move(x),
+                AsGenericExpr(ConvertTo(ix, std::move(by))), defaultRealKind);
           },
           [&](Expr<SomeReal> &&rx, BOZLiteralConstant &&by) {
-            return NumericOperation<OPR, CAN_BE_UNSIGNED>(messages,
-                std::move(x), AsGenericExpr(ConvertTo(rx, std::move(by))),
-                defaultRealKind);
+            return NumericOperation<OPR>(messages, std::move(x),
+                AsGenericExpr(ConvertTo(rx, std::move(by))), defaultRealKind);
           },
           // Error cases
           [&](Expr<SomeUnsigned> &&, auto &&) {
@@ -621,7 +613,7 @@ std::optional<Expr<SomeType>> NumericOperation(
       std::move(x.u), std::move(y.u));
 }
 
-template std::optional<Expr<SomeType>> NumericOperation<Power, false>(
+template std::optional<Expr<SomeType>> NumericOperation<Power>(
     parser::ContextualMessages &, Expr<SomeType> &&, Expr<SomeType> &&,
     int defaultRealKind);
 template std::optional<Expr<SomeType>> NumericOperation<Multiply>(
diff --git a/flang/lib/Optimizer/Builder/IntrinsicCall.cpp b/flang/lib/Optimizer/Builder/IntrinsicCall.cpp
index 70a373479837..e1c9520592de 100644
--- a/flang/lib/Optimizer/Builder/IntrinsicCall.cpp
+++ b/flang/lib/Optimizer/Builder/IntrinsicCall.cpp
@@ -138,7 +138,7 @@ static const char __ldlu_r8x2[] = "__ldlu_r8x2_";
 /// Table that drives the fir generation depending on the intrinsic or intrinsic
 /// module procedure one to one mapping with Fortran arguments. If no mapping is
 /// defined here for a generic intrinsic, genRuntimeCall will be called
-/// to look for a match in the runtime a emit a call. Note that the argument
+/// to look for a match in the runtime and emit a call. Note that the argument
 /// lowering rules for an intrinsic need to be provided only if at least one
 /// argument must not be lowered by value. In which case, the lowering rules
 /// should be provided for all the intrinsic arguments for completeness.
@@ -1062,7 +1062,7 @@ prettyPrintIntrinsicName(fir::FirOpBuilder &builder, mlir::Location loc,
                          llvm::StringRef suffix, mlir::FunctionType funcType) {
   std::string output = prefix.str();
   llvm::raw_string_ostream sstream(output);
-  if (name == "pow") {
+  if (name == "pow" || name == "pow-unsigned") {
     assert(funcType.getNumInputs() == 2 && "power operator has two arguments");
     std::string displayName{" ** "};
     sstream << mlirTypeToIntrinsicFortran(builder, funcType.getInput(0), loc,
@@ -1675,6 +1675,14 @@ static constexpr MathOperation mathOperations[] = {
      genComplexPow},
     {"pow", RTNAME_STRING(cqpowk), FuncTypeComplex16Complex16Integer8,
      genLibF128Call},
+    {"pow-unsigned", RTNAME_STRING(UPow1),
+     genFuncType<Ty::Integer<1>, Ty::Integer<1>, Ty::Integer<1>>, genLibCall},
+    {"pow-unsigned", RTNAME_STRING(UPow2),
+     genFuncType<Ty::Integer<2>, Ty::Integer<2>, Ty::Integer<2>>, genLibCall},
+    {"pow-unsigned", RTNAME_STRING(UPow4),
+     genFuncType<Ty::Integer<4>, Ty::Integer<4>, Ty::Integer<4>>, genLibCall},
+    {"pow-unsigned", RTNAME_STRING(UPow8),
+     genFuncType<Ty::Integer<8>, Ty::Integer<8>, Ty::Integer<8>>, genLibCall},
     {"remainder", "remainderf",
      genFuncType<Ty::Real<4>, Ty::Real<4>, Ty::Real<4>>, genLibCall},
     {"remainder", "remainder",
@@ -9441,6 +9449,14 @@ mlir::Value genPow(fir::FirOpBuilder &builder, mlir::Location loc,
   //       implementation and mark it 'strictfp'.
   //       Another option is to implement it in Fortran runtime library
   //       (just like matmul).
+  if (type.isUnsignedInteger()) {
+    assert(x.getType().isUnsignedInteger() && y.getType().isUnsignedInteger() &&
+           "unsigned pow requires unsigned arguments");
+    return IntrinsicLibrary{builder, loc}.genRuntimeCall("pow-unsigned", type,
+                                                         {x, y});
+  }
+  assert(!x.getType().isUnsignedInteger() && !y.getType().isUnsignedInteger() &&
+         "non-unsigned pow requires non-unsigned arguments");
   return IntrinsicLibrary{builder, loc}.genRuntimeCall("pow", type, {x, y});
 }
 
diff --git a/flang/lib/Semantics/expression.cpp b/flang/lib/Semantics/expression.cpp
index c1f01a0af94f..ea587b6831b7 100644
--- a/flang/lib/Semantics/expression.cpp
+++ b/flang/lib/Semantics/expression.cpp
@@ -3783,10 +3783,9 @@ MaybeExpr NumericBinaryHelper(
       analyzer.CheckForNullPointer();
       analyzer.CheckForAssumedRank();
       analyzer.CheckConformance();
-      constexpr bool canBeUnsigned{opr != NumericOperator::Power};
-      return NumericOperation<OPR, canBeUnsigned>(
-          context.GetContextualMessages(), analyzer.MoveExpr(0),
-          analyzer.MoveExpr(1), context.GetDefaultKind(TypeCategory::Real));
+      return NumericOperation<OPR>(context.GetContextualMessages(),
+          analyzer.MoveExpr(0), analyzer.MoveExpr(1),
+          context.GetDefaultKind(TypeCategory::Real));
     } else {
       return analyzer.TryDefinedOp(AsFortran(opr),
           "Operands of %s must be numeric; have %s and %s"_err_en_US);
diff --git a/flang/test/Lower/unsigned-ops.f90 b/flang/test/Lower/unsigned-ops.f90
index f61f10656159..13e17721ceb9 100644
--- a/flang/test/Lower/unsigned-ops.f90
+++ b/flang/test/Lower/unsigned-ops.f90
@@ -24,3 +24,29 @@ end
 !CHECK: fir.store %[[VAL_13]] to %[[VAL_2]] : !fir.ref<ui32>
 !CHECK: %[[VAL_14:.*]] = fir.load %[[VAL_2]] : !fir.ref<ui32>
 !CHECK: return %[[VAL_14]] : ui32
+
+unsigned function f02(u, v)
+  unsigned, intent(in) :: u, v
+  f02 = u ** v - 1u
+end
+
+!CHECK: func.func @_QPf02(%[[ARG0:.*]]: !fir.ref<ui32> {fir.bindc_name = "u"}, %[[ARG1:.*]]: !fir.ref<ui32> {fir.bindc_name = "v"}) -> ui32 {
+!CHECK: %[[C1_i32:.*]] = arith.constant 1 : i32
+!CHECK: %[[VAL_0:.*]] = fir.dummy_scope : !fir.dscope
+!CHECK: %[[VAL_1:.*]] = fir.alloca ui32 {bindc_name = "f02", uniq_name = "_QFf02Ef02"}
+!CHECK: %[[VAL_2:.*]] = fir.declare %[[VAL_1]] {uniq_name = "_QFf02Ef02"} : (!fir.ref<ui32>) -> !fir.ref<ui32>
+!CHECK: %[[VAL_3:.*]] = fir.declare %[[ARG0]] dummy_scope %[[VAL_0]] {fortran_attrs = #fir.var_attrs<intent_in>, uniq_name = "_QFf02Eu"} : (!fir.ref<ui32>, !fir.dscope) -> !fir.ref<ui32>
+!CHECK: %[[VAL_4:.*]] = fir.declare %[[ARG1]] dummy_scope %[[VAL_0]] {fortran_attrs = #fir.var_attrs<intent_in>, uniq_name = "_QFf02Ev"} : (!fir.ref<ui32>, !fir.dscope) -> !fir.ref<ui32>
+!CHECK: %[[VAL_5:.*]] = fir.load %[[VAL_3]] : !fir.ref<ui32>
+!CHECK: %[[VAL_6:.*]] = fir.load %[[VAL_4]] : !fir.ref<ui32>
+!CHECK: %[[VAL_7:.*]] = fir.convert %[[VAL_5]] : (ui32) -> i64
+!CHECK: %[[VAL_8:.*]] = fir.convert %[[VAL_6]] : (ui32) -> i64
+!CHECK: %[[VAL_9:.*]] = fir.call @_FortranAUPow8(%[[VAL_7]], %[[VAL_8]]) fastmath<contract> : (i64, i64) -> i64
+!CHECK: %[[VAL_10:.*]] = fir.convert %[[VAL_9]] : (i64) -> ui32
+!CHECK: %[[VAL_11:.*]] = fir.convert %[[VAL_10]] : (ui32) -> i32
+!CHECK: %[[VAL_12:.*]] = arith.subi %[[VAL_11]], %[[C1_i32]] : i32
+!CHECK: %[[VAL_13:.*]] = fir.convert %[[VAL_12]] : (i32) -> ui32
+!CHECK: fir.store %[[VAL_13]] to %[[VAL_2]] : !fir.ref<ui32>
+!CHECK: %[[VAL_14:.*]] = fir.load %[[VAL_2]] : !fir.ref<ui32>
+!CHECK: return %[[VAL_14]] : ui32
+
diff --git a/flang/test/Semantics/unsigned-errors.f90 b/flang/test/Semantics/unsigned-errors.f90
index 2e2539b40e5e..18f28f21ef8a 100644
--- a/flang/test/Semantics/unsigned-errors.f90
+++ b/flang/test/Semantics/unsigned-errors.f90
@@ -20,8 +20,7 @@ print *, 0u + 1u ! ok
 print *, 0u - 1u ! ok
 print *, 0u * 1u ! ok
 print *, 0u / 1u ! ok
-!ERROR: Operands must not be UNSIGNED
-print *, 0u ** 1u
+print *, 0u ** 1u ! ok
 
 print *, uint((0.,0.)) ! ok
 print *, uint(z'123') ! ok