diff --git a/CMakeLists.txt b/CMakeLists.txt
index e577fb38ce..0210ac01a1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -16,9 +16,6 @@ set(CMAKE_MODULE_PATH "${LAPACK_SOURCE_DIR}/CMAKE" ${CMAKE_MODULE_PATH})
 # Export all symbols on Windows when building shared libraries
 SET(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS TRUE)
 
-# Tell CMake that our Fortran sources are written in fixed format.
-set(CMAKE_Fortran_FORMAT FIXED)
-
 # Set a default build type if none was specified
 if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
   message(STATUS "Setting build type to 'Release' as none was specified.")
diff --git a/DOCS/Doxyfile b/DOCS/Doxyfile
index 43cea43b5a..92057726e8 100644
--- a/DOCS/Doxyfile
+++ b/DOCS/Doxyfile
@@ -793,6 +793,7 @@ INPUT_ENCODING         = UTF-8
 
 FILE_PATTERNS          = *.c \
                          *.f \
+                         *.f90 \
                          *.h
 
 # The RECURSIVE tag can be used to specify whether or not subdirectories should
diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt
index bd3be3c6e8..62e6ec72f1 100644
--- a/SRC/CMakeLists.txt
+++ b/SRC/CMakeLists.txt
@@ -41,6 +41,7 @@ set(ALLAUX ilaenv.f ilaenv2stage.f ieeeck.f lsamen.f iparmq.f iparam2stage.F
    ../INSTALL/slamch.f)
 
 set(SCLAUX
+   la_constants.f90
    sbdsdc.f
    sbdsqr.f sdisna.f slabad.f slacpy.f sladiv.f slae2.f  slaebz.f
    slaed0.f slaed1.f slaed2.f slaed3.f slaed4.f slaed5.f slaed6.f
@@ -59,6 +60,7 @@ set(SCLAUX
    ${SECOND_SRC})
 
 set(DZLAUX
+   la_constants.f90
    dbdsdc.f
    dbdsqr.f ddisna.f dlabad.f dlacpy.f dladiv.f dlae2.f  dlaebz.f
    dlaed0.f dlaed1.f dlaed2.f dlaed3.f dlaed4.f dlaed5.f dlaed6.f
diff --git a/SRC/Makefile b/SRC/Makefile
index 48b0a91e94..68f1495cc9 100644
--- a/SRC/Makefile
+++ b/SRC/Makefile
@@ -57,9 +57,19 @@
 TOPSRCDIR = ..
 include $(TOPSRCDIR)/make.inc
 
-.SUFFIXES: .F .o
-.F.o:
+ALLMOD = la_constants.mod
+
+.SUFFIXES: .f .F .f90 .F90 .o .mod
+%.o: %.f $(ALLMOD)
+	$(FC) $(FFLAGS) -c -o $@ $<
+%.o: %.F $(ALLMOD)
+	$(FC) $(FFLAGS) -c -o $@ $<
+%.o: %.f90 $(ALLMOD)
 	$(FC) $(FFLAGS) -c -o $@ $<
+%.o: %.F90 $(ALLMOD)
+	$(FC) $(FFLAGS) -c -o $@ $<
+.o.mod:
+	@true
 
 ALLAUX = ilaenv.o ilaenv2stage.o ieeeck.o lsamen.o xerbla.o xerbla_array.o \
    iparmq.o iparam2stage.o \
@@ -67,6 +77,7 @@ ALLAUX = ilaenv.o ilaenv2stage.o ieeeck.o lsamen.o xerbla.o xerbla_array.o \
    ../INSTALL/ilaver.o ../INSTALL/lsame.o ../INSTALL/slamch.o
 
 SCLAUX = \
+   la_constants.o \
    sbdsdc.o \
    sbdsqr.o sdisna.o slabad.o slacpy.o sladiv.o slae2.o  slaebz.o \
    slaed0.o slaed1.o slaed2.o slaed3.o slaed4.o slaed5.o slaed6.o \
@@ -85,6 +96,7 @@ SCLAUX = \
    ../INSTALL/second_$(TIMER).o
 
 DZLAUX = \
+   la_constants.o \
    dbdsdc.o \
    dbdsqr.o ddisna.o dlabad.o dlacpy.o dladiv.o dlae2.o  dlaebz.o \
    dlaed0.o dlaed1.o dlaed2.o dlaed3.o dlaed4.o dlaed5.o dlaed6.o \
@@ -571,7 +583,7 @@ FRC:
 .PHONY: clean cleanobj cleanlib
 clean: cleanobj cleanlib
 cleanobj:
-	rm -f *.o DEPRECATED/*.o
+	rm -f *.o *.mod DEPRECATED/*.o DEPRECATED/*.mod
 cleanlib:
 	rm -f $(LAPACKLIB)
 
@@ -581,3 +593,7 @@ sla_wwaddw.o: sla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
 dla_wwaddw.o: dla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
 cla_wwaddw.o: cla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
 zla_wwaddw.o: zla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
+
+# Modules
+la_constants.o: la_constants.f90
+	$(FC) $(FFLAGS) -c -o $@ $<
diff --git a/SRC/la_constants.f90 b/SRC/la_constants.f90
new file mode 100644
index 0000000000..d8146cb1bb
--- /dev/null
+++ b/SRC/la_constants.f90
@@ -0,0 +1,130 @@
+!> \brief \b LA_CONSTANTS is a module for the scaling constants for the compiled Fortran single and double precisions
+!
+!  =========== DOCUMENTATION ===========
+!
+! Online html documentation available at
+!            http://www.netlib.org/lapack/explore-html/
+!
+!  Authors:
+!  ========
+!
+!> \author Edward Anderson, Lockheed Martin
+!
+!> \date May 2016
+!
+!> \ingroup OTHERauxiliary
+!
+!> \par Contributors:
+!  ==================
+!>
+!> Weslley Pereira, University of Colorado Denver, USA
+!> Nick Papior, Technical University of Denmark, DK
+!
+!> \par Further Details:
+!  =====================
+!>
+!> \verbatim
+!>
+!>  Anderson E. (2017)
+!>  Algorithm 978: Safe Scaling in the Level 1 BLAS
+!>  ACM Trans Math Softw 44:1--28
+!>  https://doi.org/10.1145/3061665
+!>
+!>  Blue, James L. (1978)
+!>  A Portable Fortran Program to Find the Euclidean Norm of a Vector
+!>  ACM Trans Math Softw 4:15--23
+!>  https://doi.org/10.1145/355769.355771
+!>
+!> \endverbatim
+!
+module LA_CONSTANTS
+!  -- LAPACK auxiliary module (version 3.10.0) --
+!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+!     February 2021
+
+!  Standard constants for 
+   integer, parameter :: sp = kind(1.e0)
+
+   real(sp), parameter :: szero = 0.0_sp
+   real(sp), parameter :: shalf = 0.5_sp
+   real(sp), parameter :: sone = 1.0_sp
+   real(sp), parameter :: stwo = 2.0_sp
+   real(sp), parameter :: sthree = 3.0_sp
+   real(sp), parameter :: sfour = 4.0_sp
+   real(sp), parameter :: seight = 8.0_sp
+   real(sp), parameter :: sten = 10.0_sp
+   complex(sp), parameter :: czero = ( 0.0_sp, 0.0_sp )
+   complex(sp), parameter :: chalf = ( 0.5_sp, 0.0_sp )
+   complex(sp), parameter :: cone = ( 1.0_sp, 0.0_sp )
+   character*1, parameter :: sprefix = 'S'
+   character*1, parameter :: cprefix = 'C'
+
+!  Scaling constants
+   real(sp), parameter :: sulp = epsilon(0._sp)
+   real(sp), parameter :: seps = sulp * 0.5_sp
+   real(sp), parameter :: ssafmin = real(radix(0._sp),sp)**max( &
+      minexponent(0._sp)-1, &
+      1-maxexponent(0._sp) &
+   )
+   real(sp), parameter :: ssafmax = sone / ssafmin
+   real(sp), parameter :: ssmlnum = ssafmin / sulp
+   real(sp), parameter :: sbignum = ssafmax * sulp
+   real(sp), parameter :: srtmin = sqrt(ssmlnum)
+   real(sp), parameter :: srtmax = sqrt(sbignum)
+
+!  Blue's scaling constants
+   real(sp), parameter :: stsml = real(radix(0._sp), sp)**ceiling( &
+       (minexponent(0._sp) - 1) * 0.5_sp)
+   real(sp), parameter :: stbig = real(radix(0._sp), sp)**floor( &
+       (maxexponent(0._sp) - digits(0._sp) + 1) * 0.5_sp)
+!  ssml = 1/s, where s was defined in https://doi.org/10.1145/355769.355771
+   real(sp), parameter :: sssml = real(radix(0._sp), sp)**( - floor( &
+       (minexponent(0._sp) - 1) * 0.5_sp))
+!  ssml = 1/S, where S was defined in https://doi.org/10.1145/355769.355771
+   real(sp), parameter :: ssbig = real(radix(0._sp), sp)**( - ceiling( &
+       (maxexponent(0._sp) - digits(0._sp) + 1) * 0.5_sp))
+
+!  Standard constants for 
+   integer, parameter :: dp = kind(1.d0)
+
+   real(dp), parameter :: dzero = 0.0_dp
+   real(dp), parameter :: dhalf = 0.5_dp
+   real(dp), parameter :: done = 1.0_dp
+   real(dp), parameter :: dtwo = 2.0_dp
+   real(dp), parameter :: dthree = 3.0_dp
+   real(dp), parameter :: dfour = 4.0_dp
+   real(dp), parameter :: deight = 8.0_dp
+   real(dp), parameter :: dten = 10.0_dp
+   complex(dp), parameter :: zzero = ( 0.0_dp, 0.0_dp )
+   complex(dp), parameter :: zhalf = ( 0.5_dp, 0.0_dp )
+   complex(dp), parameter :: zone = ( 1.0_dp, 0.0_dp )
+   character*1, parameter :: dprefix = 'D'
+   character*1, parameter :: zprefix = 'Z'
+
+!  Scaling constants
+   real(dp), parameter :: dulp = epsilon(0._dp)
+   real(dp), parameter :: deps = dulp * 0.5_dp
+   real(dp), parameter :: dsafmin = real(radix(0._dp),dp)**max( &
+      minexponent(0._dp)-1, &
+      1-maxexponent(0._dp) &
+   )
+   real(dp), parameter :: dsafmax = done / dsafmin
+   real(dp), parameter :: dsmlnum = dsafmin / dulp
+   real(dp), parameter :: dbignum = dsafmax * dulp
+   real(dp), parameter :: drtmin = sqrt(dsmlnum)
+   real(dp), parameter :: drtmax = sqrt(dbignum)
+
+!  Blue's scaling constants
+   real(dp), parameter :: dtsml = real(radix(0._dp), dp)**ceiling( &
+       (minexponent(0._dp) - 1) * 0.5_dp)
+   real(dp), parameter :: dtbig = real(radix(0._dp), dp)**floor( &
+       (maxexponent(0._dp) - digits(0._dp) + 1) * 0.5_dp)
+!  ssml = 1/s, where s was defined in https://doi.org/10.1145/355769.355771
+   real(dp), parameter :: dssml = real(radix(0._dp), dp)**( - floor( &
+       (minexponent(0._dp) - 1) * 0.5_dp))
+!  ssml = 1/S, where S was defined in https://doi.org/10.1145/355769.355771
+   real(dp), parameter :: dsbig = real(radix(0._dp), dp)**( - ceiling( &
+       (maxexponent(0._dp) - digits(0._dp) + 1) * 0.5_dp))
+
+end module LA_CONSTANTS