diff --git a/flang/include/flang/Evaluate/real.h b/flang/include/flang/Evaluate/real.h
--- a/flang/include/flang/Evaluate/real.h
+++ b/flang/include/flang/Evaluate/real.h
@@ -119,14 +119,21 @@
       const Real &, Rounding rounding = defaultRounding) const;
 
   ValueWithRealFlags<Real> SQRT(Rounding rounding = defaultRounding) const;
-
   // NEAREST(), IEEE_NEXT_AFTER(), IEEE_NEXT_UP(), and IEEE_NEXT_DOWN()
   ValueWithRealFlags<Real> NEAREST(bool upward) const;
-
   // HYPOT(x,y)=SQRT(x**2 + y**2) computed so as to avoid spurious
   // intermediate overflows.
   ValueWithRealFlags<Real> HYPOT(
       const Real &, Rounding rounding = defaultRounding) const;
+  // DIM(X,Y) = MAX(X-Y, 0)
+  ValueWithRealFlags<Real> DIM(
+      const Real &, Rounding rounding = defaultRounding) const;
+  // MOD(x,y) = x - AINT(x/y)*y
+  // MODULO(x,y) = x - FLOOR(x/y)*y
+  ValueWithRealFlags<Real> MOD(
+      const Real &, Rounding rounding = defaultRounding) const;
+  ValueWithRealFlags<Real> MODULO(
+      const Real &, Rounding rounding = defaultRounding) const;
 
   // DIM(X,Y) = MAX(X-Y, 0)
   ValueWithRealFlags<Real> DIM(
diff --git a/flang/lib/Evaluate/fold-real.cpp b/flang/lib/Evaluate/fold-real.cpp
--- a/flang/lib/Evaluate/fold-real.cpp
+++ b/flang/lib/Evaluate/fold-real.cpp
@@ -72,7 +72,7 @@
   } else if (name == "amax0" || name == "amin0" || name == "amin1" ||
       name == "amax1" || name == "dmin1" || name == "dmax1") {
     return RewriteSpecificMINorMAX(context, std::move(funcRef));
-  } else if (name == "atan" || name == "atan2" || name == "mod") {
+  } else if (name == "atan" || name == "atan2") {
     std::string localName{name == "atan" ? "atan2" : name};
     CHECK(args.size() == 2);
     if (auto callable{GetHostRuntimeWrapper<T, T, T>(localName)}) {
@@ -159,6 +159,35 @@
         RelationalOperator::GT, T::Scalar::HUGE().Negate());
   } else if (name == "merge") {
     return FoldMerge<T>(context, std::move(funcRef));
+  } else if (name == "min") {
+    return FoldMINorMAX(context, std::move(funcRef), Ordering::Less);
+  } else if (name == "minval") {
+    return FoldMaxvalMinval<T>(
+        context, std::move(funcRef), RelationalOperator::LT, T::Scalar::HUGE());
+  } else if (name == "mod") {
+    CHECK(args.size() == 2);
+    return FoldElementalIntrinsic<T, T, T>(context, std::move(funcRef),
+        ScalarFunc<T, T, T>(
+            [&context](const Scalar<T> &x, const Scalar<T> &y) -> Scalar<T> {
+              auto result{x.MOD(y)};
+              if (result.flags.test(RealFlag::DivideByZero)) {
+                context.messages().Say(
+                    "second argument to MOD must not be zero"_warn_en_US);
+              }
+              return result.value;
+            }));
+  } else if (name == "modulo") {
+    CHECK(args.size() == 2);
+    return FoldElementalIntrinsic<T, T, T>(context, std::move(funcRef),
+        ScalarFunc<T, T, T>(
+            [&context](const Scalar<T> &x, const Scalar<T> &y) -> Scalar<T> {
+              auto result{x.MODULO(y)};
+              if (result.flags.test(RealFlag::DivideByZero)) {
+                context.messages().Say(
+                    "second argument to MODULO must not be zero"_warn_en_US);
+              }
+              return result.value;
+            }));
   } else if (name == "nearest") {
     if (const auto *sExpr{UnwrapExpr<Expr<SomeReal>>(args[1])}) {
       return common::visit(
@@ -184,11 +213,6 @@
           },
           sExpr->u);
     }
-  } else if (name == "min") {
-    return FoldMINorMAX(context, std::move(funcRef), Ordering::Less);
-  } else if (name == "minval") {
-    return FoldMaxvalMinval<T>(
-        context, std::move(funcRef), RelationalOperator::LT, T::Scalar::HUGE());
   } else if (name == "product") {
     auto one{Scalar<T>::FromInteger(value::Integer<8>{1}).value};
     return FoldProduct<T>(context, std::move(funcRef), one);
diff --git a/flang/lib/Evaluate/intrinsics-library.cpp b/flang/lib/Evaluate/intrinsics-library.cpp
--- a/flang/lib/Evaluate/intrinsics-library.cpp
+++ b/flang/lib/Evaluate/intrinsics-library.cpp
@@ -230,7 +230,6 @@
       FolderFactory<F, F{std::log}>::Create("log"),
       FolderFactory<F, F{std::log10}>::Create("log10"),
       FolderFactory<F, F{std::lgamma}>::Create("log_gamma"),
-      FolderFactory<F2, F2{std::fmod}>::Create("mod"),
       FolderFactory<F2, F2{std::pow}>::Create("pow"),
       FolderFactory<F, F{std::sin}>::Create("sin"),
       FolderFactory<F, F{std::sinh}>::Create("sinh"),
diff --git a/flang/lib/Evaluate/real.cpp b/flang/lib/Evaluate/real.cpp
--- a/flang/lib/Evaluate/real.cpp
+++ b/flang/lib/Evaluate/real.cpp
@@ -422,6 +422,32 @@
   return result;
 }
 
+// MOD(x,y) = x - AINT(x/y)*y
+template <typename W, int P>
+ValueWithRealFlags<Real<W, P>> Real<W, P>::MOD(
+    const Real &y, Rounding rounding) const {
+  ValueWithRealFlags<Real> result;
+  Real quotient{Divide(y, rounding).AccumulateFlags(result.flags)};
+  Real toInt{quotient.ToWholeNumber(common::RoundingMode::ToZero)
+                 .AccumulateFlags(result.flags)};
+  Real product{toInt.Multiply(y, rounding).AccumulateFlags(result.flags)};
+  result.value = Subtract(product, rounding).AccumulateFlags(result.flags);
+  return result;
+}
+
+// MODULO(x,y) = x - FLOOR(x/y)*y
+template <typename W, int P>
+ValueWithRealFlags<Real<W, P>> Real<W, P>::MODULO(
+    const Real &y, Rounding rounding) const {
+  ValueWithRealFlags<Real> result;
+  Real quotient{Divide(y, rounding).AccumulateFlags(result.flags)};
+  Real toInt{quotient.ToWholeNumber(common::RoundingMode::Down)
+                 .AccumulateFlags(result.flags)};
+  Real product{toInt.Multiply(y, rounding).AccumulateFlags(result.flags)};
+  result.value = Subtract(product, rounding).AccumulateFlags(result.flags);
+  return result;
+}
+
 template <typename W, int P>
 ValueWithRealFlags<Real<W, P>> Real<W, P>::DIM(
     const Real &y, Rounding rounding) const {
diff --git a/flang/test/Evaluate/fold-mod.f90 b/flang/test/Evaluate/fold-mod.f90
new file mode 100644
--- /dev/null
+++ b/flang/test/Evaluate/fold-mod.f90
@@ -0,0 +1,24 @@
+! RUN: %python %S/test_folding.py %s %flang_fc1
+! Tests folding of integer and real MOD and MODULO
+module m1
+  logical, parameter :: test_mod_i1 = mod(8, 5) == 3
+  logical, parameter :: test_mod_i2 = mod(-8, 5) == -3
+  logical, parameter :: test_mod_i3 = mod(8, -5) == 3
+  logical, parameter :: test_mod_i4 = mod(-8, -5) == -3
+
+  logical, parameter :: test_mod_r1 = mod(3., 2.) == 1.
+  logical, parameter :: test_mod_r2 = mod(8., 5.) == 3.
+  logical, parameter :: test_mod_r3 = mod(-8., 5.) == -3.
+  logical, parameter :: test_mod_r4 = mod(8., -5.) == 3.
+  logical, parameter :: test_mod_r5 = mod(-8., -5.) == -3.
+
+  logical, parameter :: test_modulo_i1 = modulo(8, 5) == 3
+  logical, parameter :: test_modulo_i2 = modulo(-8, 5) == 2
+  logical, parameter :: test_modulo_i3 = modulo(8, -5) == -2
+  logical, parameter :: test_modulo_i4 = modulo(-8, -5) == -3
+
+  logical, parameter :: test_modulo_r1 = modulo(8., 5.) == 3.
+  logical, parameter :: test_modulo_r2 = modulo(-8., 5.) == 2.
+  logical, parameter :: test_modulo_r3 = modulo(8., -5.) == -2.
+  logical, parameter :: test_modulo_r4 = modulo(-8., -5.) == -3.
+end module
diff --git a/flang/test/Evaluate/folding04.f90 b/flang/test/Evaluate/folding04.f90
--- a/flang/test/Evaluate/folding04.f90
+++ b/flang/test/Evaluate/folding04.f90
@@ -32,6 +32,9 @@
   !WARN: warning: invalid argument on intrinsic function
   real(4), parameter :: nan_r4_acos5 = acos(r4_pinf)
   TEST_ISNAN(nan_r4_acos5)
+  !WARN: warning: second argument to MOD must not be zero
+  real(4), parameter :: nan_r4_mod = mod(3.5, 0.)
+  TEST_ISNAN(nan_r4_mod)
 
   !WARN: warning: overflow on intrinsic function
   logical, parameter :: test_exp_overflow = exp(256._4).EQ.r4_pinf