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
@@ -280,15 +280,34 @@
     // SQRT(+Inf) == +Inf
     result.value = Infinity(false);
   } else {
-    // Slow but reliable bit-at-a-time method.  Start with a clear significand
-    // and half the unbiased exponent, and then try to set significand bits
-    // in descending order of magnitude without exceeding the exact result.
     int expo{UnbiasedExponent()};
-    if (IsSubnormal()) {
-      expo -= GetFraction().LEADZ();
+    if (expo < -1 || expo > 1) {
+      // Reduce the range to [0.5 .. 4.0) by dividing by an integral power
+      // of four to avoid trouble with very large and very small values
+      // (esp. truncation of subnormals).
+      // SQRT(2**(2a) * x) = SQRT(2**(2a)) * SQRT(x) = 2**a * SQRT(x)
+      Real scaled;
+      int adjust{expo / 2};
+      scaled.Normalize(false, expo - 2 * adjust + exponentBias, GetFraction());
+      result = scaled.SQRT(rounding);
+      result.value.Normalize(false,
+          result.value.UnbiasedExponent() + adjust + exponentBias,
+          result.value.GetFraction());
+      return result;
     }
+    // Compute the square root of the reduced value with the slow but
+    // reliable bit-at-a-time method.  Start with a clear significand and
+    // half of the unbiased exponent, and then try to set significand bits
+    // in descending order of magnitude without exceeding the exact result.
     expo = expo / 2 + exponentBias;
     result.value.Normalize(false, expo, Fraction::MASKL(1));
+    Real initialSq{result.value.Multiply(result.value).value};
+    if (Compare(initialSq) == Relation::Less) {
+      // Initial estimate is too large; this can happen for values just
+      // under 1.0.
+      --expo;
+      result.value.Normalize(false, expo, Fraction::MASKL(1));
+    }
     for (int bit{significandBits - 1}; bit >= 0; --bit) {
       Word word{result.value.word_};
       result.value.word_ = word.IBSET(bit);
@@ -299,10 +318,10 @@
         result.value.word_ = word;
       }
     }
-    // The computed square root, when squared, has a square that's not greater
-    // than the original argument.  Check this square against the square of the
-    // next Real value, and return that one if its square is closer in magnitude
-    // to the original argument.
+    // The computed square root has a square that's not greater than the
+    // original argument.  Check this square against the square of the next
+    // larger Real and return that one if its square is closer in magnitude to
+    // the original argument.
     Real resultSq{result.value.Multiply(result.value).value};
     Real diff{Subtract(resultSq).value.ABS()};
     if (diff.IsZero()) {
diff --git a/flang/test/Evaluate/folding28.f90 b/flang/test/Evaluate/folding28.f90
--- a/flang/test/Evaluate/folding28.f90
+++ b/flang/test/Evaluate/folding28.f90
@@ -27,7 +27,7 @@
   logical, parameter :: test_sqr_sqrt_t8 = sqr_sqrt_t8 == t8
   ! max subnormal
   real(8), parameter :: maxs8 = z'000fffffffffffff'
-  real(8), parameter :: sqrt_maxs8 = sqrt(maxs8), sqrt_maxs8z = z'2000000000000000'
+  real(8), parameter :: sqrt_maxs8 = sqrt(maxs8), sqrt_maxs8z = z'1fffffffffffffff'
   logical, parameter :: test_sqrt_maxs8 = sqrt_maxs8 == sqrt_maxs8z
   ! min subnormal
   real(8), parameter :: mins8 = z'1'
@@ -35,5 +35,13 @@
   logical, parameter :: test_sqrt_mins8 = sqrt_mins8 == sqrt_mins8z
   real(8), parameter :: sqr_sqrt_mins8 = sqrt_mins8 * sqrt_mins8
   logical, parameter :: test_sqr_sqrt_mins8 = sqr_sqrt_mins8 == mins8
+  ! regression tests: cases near 1.
+  real(4), parameter :: sqrt_under1 = sqrt(.96875)
+  logical, parameter :: test_sqrt_under1 = sqrt_under1 == .984250962734222412109375
+  ! oddball case: the value before 1. is also its own sqrt, but not its own square
+  real(4), parameter :: before_1 = z'3f7fffff' ! .999999940395355224609375
+  real(4), parameter :: sqrt_before_1 = sqrt(before_1)
+  logical, parameter :: test_before_1 = sqrt_before_1 == before_1
+  real(4), parameter :: sq_sqrt_before_1 = sqrt_before_1 * sqrt_before_1
+  logical, parameter :: test_sq_before_1 = sq_sqrt_before_1 < before_1
 end module
-