Index: libgfortran/config/fpu-387.h
===================================================================
--- libgfortran/config/fpu-387.h	(revision 212239)
+++ libgfortran/config/fpu-387.h	(working copy)
@@ -64,6 +64,11 @@ has_sse (void)
 
 #define _FPU_RC_MASK    0x3
 
+/* Enable flush to zero mode.  */
+
+#define MXCSR_FTZ (1 << 15)
+
+
 /* This structure corresponds to the layout of the block
    written by FSTENV.  */
 typedef struct
@@ -458,3 +463,44 @@ set_fpu_state (void *state)
     __asm__ __volatile__ ("%vldmxcsr\t%0" : : "m" (envp->__mxcsr));
 }
 
+
+int
+support_fpu_underflow_control (int kind)
+{
+  return (has_sse() && (kind == 4 || kind == 8)) ? 1 : 0;
+}
+
+
+int
+get_fpu_underflow_mode (void)
+{
+  unsigned int cw_sse;
+
+  if (!has_sse())
+    return 1;
+
+  __asm__ __volatile__ ("%vstmxcsr\t%0" : "=m" (cw_sse));
+
+  /* Return 0 for abrupt underflow (flush to zero), 1 for gradual underflow.  */
+  return (cw_sse & MXCSR_FTZ) ? 0 : 1;
+}
+
+
+void
+set_fpu_underflow_mode (int gradual)
+{
+  unsigned int cw_sse;
+
+  if (!has_sse())
+    return;
+
+  __asm__ __volatile__ ("%vstmxcsr\t%0" : "=m" (cw_sse));
+
+  if (gradual)
+    cw_sse &= ~MXCSR_FTZ;
+  else
+    cw_sse |= MXCSR_FTZ;
+
+  __asm__ __volatile__ ("%vldmxcsr\t%0" : : "m" (cw_sse));
+}
+
Index: libgfortran/config/fpu-aix.h
===================================================================
--- libgfortran/config/fpu-aix.h	(revision 212239)
+++ libgfortran/config/fpu-aix.h	(working copy)
@@ -418,3 +418,23 @@ set_fpu_state (void *state)
   fesetenv (state);
 }
 
+
+int
+support_fpu_underflow_control (int kind)
+{
+  return 0;
+}
+
+
+int
+get_fpu_underflow_mode (void)
+{
+  return 0;
+}
+
+
+void
+set_fpu_underflow_mode (int gradual __attribute__((unused)))
+{
+}
+
Index: libgfortran/config/fpu-generic.h
===================================================================
--- libgfortran/config/fpu-generic.h	(revision 212239)
+++ libgfortran/config/fpu-generic.h	(working copy)
@@ -75,3 +75,24 @@ void
 set_fpu_rounding_mode (int round __attribute__((unused)))
 {
 }
+
+
+int
+support_fpu_underflow_control (int kind)
+{
+  return 0;
+}
+
+
+int
+get_fpu_underflow_mode (void)
+{
+  return 0;
+}
+
+
+void
+set_fpu_underflow_mode (int gradual __attribute__((unused)))
+{
+}
+
Index: libgfortran/config/fpu-glibc.h
===================================================================
--- libgfortran/config/fpu-glibc.h	(revision 212239)
+++ libgfortran/config/fpu-glibc.h	(working copy)
@@ -432,3 +432,23 @@ set_fpu_state (void *state)
   fesetenv (state);
 }
 
+
+int
+support_fpu_underflow_control (int kind)
+{
+  return 0;
+}
+
+
+int
+get_fpu_underflow_mode (void)
+{
+  return 0;
+}
+
+
+void
+set_fpu_underflow_mode (int gradual __attribute__((unused)))
+{
+}
+
Index: libgfortran/config/fpu-sysv.h
===================================================================
--- libgfortran/config/fpu-sysv.h	(revision 212239)
+++ libgfortran/config/fpu-sysv.h	(working copy)
@@ -468,3 +468,23 @@ set_fpu_state (void *s)
   fpsetround (state->round);
 }
 
+
+int
+support_fpu_underflow_control (int kind)
+{
+  return 0;
+}
+
+
+int
+get_fpu_underflow_mode (void)
+{
+  return 0;
+}
+
+
+void
+set_fpu_underflow_mode (int gradual __attribute__((unused)))
+{
+}
+
Index: libgfortran/ieee/ieee_arithmetic.F90
===================================================================
--- libgfortran/ieee/ieee_arithmetic.F90	(revision 212239)
+++ libgfortran/ieee/ieee_arithmetic.F90	(working copy)
@@ -349,6 +349,29 @@ module IEEE_ARITHMETIC
     end function
   end interface
 
+  ! IEEE_SUPPORT_UNDERFLOW_CONTROL
+
+  interface IEEE_SUPPORT_UNDERFLOW_CONTROL
+    module procedure IEEE_SUPPORT_UNDERFLOW_CONTROL_4, &
+                     IEEE_SUPPORT_UNDERFLOW_CONTROL_8, &
+#ifdef HAVE_GFC_REAL_10
+                     IEEE_SUPPORT_UNDERFLOW_CONTROL_10, &
+#endif
+#ifdef HAVE_GFC_REAL_16
+                     IEEE_SUPPORT_UNDERFLOW_CONTROL_16, &
+#endif
+                     IEEE_SUPPORT_UNDERFLOW_CONTROL_NOARG
+  end interface
+  public :: IEEE_SUPPORT_UNDERFLOW_CONTROL
+  
+  ! Interface to the FPU-specific function
+  interface
+    pure integer function support_underflow_control_helper(kind) &
+        bind(c, name="_gfortrani_support_fpu_underflow_control")
+      integer, intent(in), value :: kind
+    end function
+  end interface
+
 ! IEEE_SUPPORT_* generic functions
 
 #if defined(HAVE_GFC_REAL_10) && defined(HAVE_GFC_REAL_16)
@@ -373,7 +396,6 @@ SUPPORTGENERIC(IEEE_SUPPORT_IO)
 SUPPORTGENERIC(IEEE_SUPPORT_NAN)
 SUPPORTGENERIC(IEEE_SUPPORT_SQRT)
 SUPPORTGENERIC(IEEE_SUPPORT_STANDARD)
-SUPPORTGENERIC(IEEE_SUPPORT_UNDERFLOW_CONTROL)
 
 contains
 
@@ -560,7 +582,6 @@ contains
   subroutine IEEE_GET_ROUNDING_MODE (ROUND_VALUE)
     implicit none
     type(IEEE_ROUND_TYPE), intent(out) :: ROUND_VALUE
-    integer :: i
 
     interface
       integer function helper() &
@@ -568,9 +589,7 @@ contains
       end function
     end interface
 
-    ! FIXME: Use intermediate variable i to avoid triggering PR59023
-    i = helper()
-    ROUND_VALUE = IEEE_ROUND_TYPE(i)
+    ROUND_VALUE = IEEE_ROUND_TYPE(helper())
   end subroutine
 
 
@@ -596,10 +615,14 @@ contains
   subroutine IEEE_GET_UNDERFLOW_MODE (GRADUAL)
     implicit none
     logical, intent(out) :: GRADUAL
-    ! We do not support getting/setting underflow mode yet. We still
-    ! provide the procedures to avoid link-time error if a user program
-    ! uses it protected by a call to IEEE_SUPPORT_UNDERFLOW_CONTROL
-    call abort
+
+    interface
+      integer function helper() &
+        bind(c, name="_gfortrani_get_fpu_underflow_mode")
+      end function
+    end interface
+
+    GRADUAL = (helper() /= 0)
   end subroutine
 
 
@@ -608,10 +631,15 @@ contains
   subroutine IEEE_SET_UNDERFLOW_MODE (GRADUAL)
     implicit none
     logical, intent(in) :: GRADUAL
-    ! We do not support getting/setting underflow mode yet. We still
-    ! provide the procedures to avoid link-time error if a user program
-    ! uses it protected by a call to IEEE_SUPPORT_UNDERFLOW_CONTROL
-    call abort
+
+    interface
+      subroutine helper(val) &
+          bind(c, name="_gfortrani_set_fpu_underflow_mode")
+        integer, value :: val
+      end subroutine
+    end interface
+
+    call helper(merge(1, 0, GRADUAL))
   end subroutine
 
 ! IEEE_SUPPORT_ROUNDING
@@ -658,6 +686,46 @@ contains
 #endif
   end function
 
+! IEEE_SUPPORT_UNDERFLOW_CONTROL
+
+  pure logical function IEEE_SUPPORT_UNDERFLOW_CONTROL_4 (X) result(res)
+    implicit none
+    real(kind=4), intent(in) :: X
+    res = (support_underflow_control_helper(4) /= 0)
+  end function
+
+  pure logical function IEEE_SUPPORT_UNDERFLOW_CONTROL_8 (X) result(res)
+    implicit none
+    real(kind=8), intent(in) :: X
+    res = (support_underflow_control_helper(8) /= 0)
+  end function
+
+#ifdef HAVE_GFC_REAL_10
+  pure logical function IEEE_SUPPORT_UNDERFLOW_CONTROL_10 (X) result(res)
+    implicit none
+    real(kind=10), intent(in) :: X
+    res = .false.
+  end function
+#endif
+
+#ifdef HAVE_GFC_REAL_16
+  pure logical function IEEE_SUPPORT_UNDERFLOW_CONTROL_16 (X) result(res)
+    implicit none
+    real(kind=16), intent(in) :: X
+    res = .false.
+  end function
+#endif
+
+  pure logical function IEEE_SUPPORT_UNDERFLOW_CONTROL_NOARG () result(res)
+    implicit none
+#if defined(HAVE_GFC_REAL_10) || defined(HAVE_GFC_REAL_16)
+    res = .false.
+#else
+    res = (support_underflow_control_helper(4) /= 0 &
+           .and. support_underflow_control_helper(8) /= 0)
+#endif
+  end function
+
 ! IEEE_SUPPORT_* functions
 
 #define SUPPORTMACRO(NAME, INTKIND, VALUE) \
@@ -801,17 +869,4 @@ SUPPORTMACRO_NOARG(IEEE_SUPPORT_STANDARD
 SUPPORTMACRO_NOARG(IEEE_SUPPORT_STANDARD,.true.)
 #endif
 
-! IEEE_SUPPORT_UNDERFLOW_CONTROL
-
-SUPPORTMACRO(IEEE_SUPPORT_UNDERFLOW_CONTROL,4,.false.)
-SUPPORTMACRO(IEEE_SUPPORT_UNDERFLOW_CONTROL,8,.false.)
-#ifdef HAVE_GFC_REAL_10
-SUPPORTMACRO(IEEE_SUPPORT_UNDERFLOW_CONTROL,10,.false.)
-#endif
-#ifdef HAVE_GFC_REAL_16
-SUPPORTMACRO(IEEE_SUPPORT_UNDERFLOW_CONTROL,16,.false.)
-#endif
-SUPPORTMACRO_NOARG(IEEE_SUPPORT_UNDERFLOW_CONTROL,.false.)
-
-
 end module IEEE_ARITHMETIC
Index: libgfortran/libgfortran.h
===================================================================
--- libgfortran/libgfortran.h	(revision 212239)
+++ libgfortran/libgfortran.h	(working copy)
@@ -787,6 +787,15 @@ internal_proto(get_fpu_state);
 extern void set_fpu_state (void *);
 internal_proto(set_fpu_state);
 
+extern int get_fpu_underflow_mode (void);
+internal_proto(get_fpu_underflow_mode);
+
+extern void set_fpu_underflow_mode (int);
+internal_proto(set_fpu_underflow_mode);
+
+extern int support_fpu_underflow_control (int);
+internal_proto(support_fpu_underflow_control);
+
 /* memory.c */
 
 extern void *xmalloc (size_t) __attribute__ ((malloc));
Index: gcc/testsuite/gfortran.dg/ieee/underflow_1.f90
===================================================================
--- gcc/testsuite/gfortran.dg/ieee/underflow_1.f90	(revision 0)
+++ gcc/testsuite/gfortran.dg/ieee/underflow_1.f90	(working copy)
@@ -0,0 +1,75 @@
+! { dg-do run }
+! { dg-additional-options "-O0" }
+! { dg-additional-options "-msse -mfpmath=sse" { target { i?86-*-* x86_64-*-* } } }
+
+program test_underflow_control
+  use ieee_arithmetic
+  use iso_fortran_env
+
+  interface use_real
+    procedure use_real_4, use_real_8
+  end interface use_real
+
+  logical l
+  real :: x
+  double precision :: y
+  integer, parameter :: kx = kind(x), ky = kind(y)
+
+  if (ieee_support_underflow_control(x)) then
+
+    x = tiny(x)
+    call ieee_set_underflow_mode(.true.)
+    call use_real(x)
+    x = x / 2000._kx
+    call use_real(x)
+    if (x <= 0) call abort
+    call ieee_get_underflow_mode(l)
+    if (.not. l) call abort
+
+    x = tiny(x)
+    call ieee_set_underflow_mode(.false.)
+    call use_real(x)
+    x = x / 2000._kx
+    call use_real(x)
+    if (x > 0) call abort
+    call ieee_get_underflow_mode(l)
+    if (l) call abort
+
+  end if
+
+  if (ieee_support_underflow_control(y)) then
+
+    y = tiny(y)
+    call ieee_set_underflow_mode(.true.)
+    call use_real(y)
+    y = y / 2000._ky
+    call use_real(y)
+    if (y <= 0) call abort
+    call ieee_get_underflow_mode(l)
+    if (.not. l) call abort
+
+    y = tiny(y)
+    call ieee_set_underflow_mode(.false.)
+    call use_real(y)
+    y = y / 2000._ky
+    call use_real(y)
+    if (y > 0) call abort
+    call ieee_get_underflow_mode(l)
+    if (l) call abort
+
+  end if
+
+contains
+
+  ! Interface to a routine that avoids calculations to be optimized out,
+  ! making it appear that we use the result
+  subroutine use_real_4(x)
+    real :: x
+    if (x == 123456.789) print *, "toto"
+  end subroutine
+  subroutine use_real_8(x)
+    double precision :: x
+    if (x == 123456.789) print *, "toto"
+  end subroutine
+
+end program