diff --git a/README.md b/README.md
index 1f43cf2..74018a8 100644
--- a/README.md
+++ b/README.md
@@ -41,11 +41,11 @@ produces output that includes actual program syntax:
 ```fortran
     f = (x**2)/2 ! <-- scalar function
     v = x        ! <-- vector function
-.SSS. (v .dot. .grad. f) * dV =  .3333333330205934
-.SSS. (     f * .div. v) * dV =  .16666666739857125
-   -.SS. (f .x. (v .dot. dA)) = -.5000000004191649
+.SSS. (v .dot. .grad. f) * dV =  .33333322
+.SSS. (     f * .div. v) * dV =  .16666671
+   -.SS. (f .x. (v .dot. dA)) = -.4999999
 ----------------------------------------------------
-                          sum = -.2220446049250313E-15 (residual)
+                          sum = .29802323E-7 (residual)
 ```
 where the small residual evidences a highly accurate approximation.
 
diff --git a/doc/uml/api.md b/doc/uml/api.md
index cf954a1..fc98007 100644
--- a/doc/uml/api.md
+++ b/doc/uml/api.md
@@ -12,16 +12,16 @@ vector_3D_t <|-- gradient_3D_t : is a
 class scalar_3D_t{
     operator(.grad.) (scalar_3D_t) gradient_3D_t
     operator(+) (scalar_3D_t, scalar_3D_t) scalar_3D_t
-    operator(*) (double precision, scalar_3D_t) scalar_3D_t
-    operator(*) (scalar_3D_t, double precision) scalar_3D_t
+    operator(*) (real, scalar_3D_t) scalar_3D_t
+    operator(*) (scalar_3D_t, real) scalar_3D_t
     operator(*) (integer, scalar_3D_t) scalar_3D_t
     operator(*) (scalar_3D_t, integer) scalar_3D_t
-    values() double precision(:,:,:)
-    grid(component : integer, coordinate : integer) double precision(:)
+    values() real(:,:,:)
+    grid(component : integer, coordinate : integer) real(:)
     consistent() : logical 
-    to_faces(direction : integer) double precision(:,:,:)
+    to_faces(direction : integer) real(:,:,:)
     to_file(name : character(len=:)) file_t
-    scalar_3D_t(initializer : scalar_3D_initializer_i, order : integer, cells : integer(:), x_min : double precision(:), x_max : double precision(:)) scalar_3D_t
+    scalar_3D_t(initializer : scalar_3D_initializer_i, order : integer, cells : integer(:), x_min : real(:), x_max : real(:)) scalar_3D_t
     scalar_3D_t(initializer : scalar_3D_initializer_i, mold : scalar_3D_t) scalar_3D_t
 }
 
@@ -30,18 +30,18 @@ class vector_3D_t{
     operator(.dot.) (vector_3D_t, vector_3D_t) : scalar_3D_t
     operator(*)(scalar_3D_t, vector_3D_t) vector_3D_t
     operator(*)(vector_3D_t, scalar_3D_t) vector_3D_t
-    grid(component : integer, coordinate : integer) double precision(:)
-    values(direction : integer) double precision(:,:,:)
+    grid(component : integer, coordinate : integer) real(:)
+    values(direction : integer) real(:,:,:)
     consistent() : logical 
-    to_centers_extended() double precision(:,:,:,:)
+    to_centers_extended() real(:,:,:,:)
     to_file(name : character(len=:)) file_t
-    vector_3D_t(initializer : vector_3D_initializer_i, order : integer, cells : integer(:), x_min : double precision(:), x_max : double precision(:))
+    vector_3D_t(initializer : vector_3D_initializer_i, order : integer, cells : integer(:), x_min : real(:), x_max : real(:))
     vector_3D_t(initializer : vector_3D_initializer_i, mold : vector_3D_t)
     vector_3D_t(initializer : vector_3D_initializer_i, mold : scalar_3D_t)
 }
 
 class divergence_3D_t{
     assignment(=)(scalar_3D_t, divergence_3D_t)
-    divergence_3D_t(initializer : vector_3D_initializer_i, order : integer, cells : integer(:), x_min : double precision(:), x_max : double precision(:))
+    divergence_3D_t(initializer : vector_3D_initializer_i, order : integer, cells : integer(:), x_min : real(:), x_max : real(:))
     divergence_3D_t(initializer : vector_3D_initializer_i, mold : vector_3D_t)
 }
diff --git a/doc/uml/class-diagram.md b/doc/uml/class-diagram.md
index 4fba59f..da90880 100644
--- a/doc/uml/class-diagram.md
+++ b/doc/uml/class-diagram.md
@@ -45,21 +45,21 @@ class vector_1D_t{
 }
 
 class gradient_1D_t{
-   - weights : double precision[]
+   - weights : real[]
 }
 
 class mimetic_matrix_1D_t{
- - upper_ :: double precision[]
- - inner_ :: double precision[]
- - lower_ :: double precision[]
+ - upper_ :: real[]
+ - inner_ :: real[]
+ - lower_ :: real[]
 }
 
 class gradient_operator_1D_t{
- + operator(.x.) double precision[]
- + assemble() double precision[] "2D array"
+ + operator(.x.) real[]
+ + assemble() real[] "2D array"
 }
 
 class divergence_operator_1D_t{
- + operator(.x.) double precision[]
- + assemble() double precision[] "2D array"
+ + operator(.x.) real[]
+ + assemble() real[] "2D array"
 }
diff --git a/example/2D-advection-diffusion.F90 b/example/2D-advection-diffusion.F90
index 6a3581d..69a474c 100644
--- a/example/2D-advection-diffusion.F90
+++ b/example/2D-advection-diffusion.F90
@@ -9,10 +9,10 @@ module fields_m
 contains
 
   pure function scalar_field(x,y) result(gaussian)
-    double precision, intent(in) :: x(:), y(:)
-    double precision gaussian(size(x),size(y))
-    double precision, parameter :: pi = acos(-1D0)
-    double precision, parameter :: x0 = -pi/2, y0 = -pi/2, sigma = pi/8
+    real, intent(in) :: x(:), y(:)
+    real gaussian(size(x),size(y))
+    real, parameter :: pi = acos(-1E0)
+    real, parameter :: x0 = -pi/2, y0 = -pi/2, sigma = pi/8
     do concurrent(integer :: j=1:size(y)) default(none) shared(x,y,gaussian)
       associate(r => sqrt((x-x0)**2 + (y(j)-y0)**2))
         gaussian(:,j) = exp(-(r**2)/(2*sigma**2))
@@ -21,8 +21,8 @@ pure function scalar_field(x,y) result(gaussian)
   end function
 
   pure function taylor_green_velocity(x,y) result(velocity)
-    double precision, intent(in) :: x(:), y(:)
-    double precision velocity(size(x),size(y),space_dimension)
+    real, intent(in) :: x(:), y(:)
+    real velocity(size(x),size(y),space_dimension)
     do concurrent(integer :: i=1:size(x), j=1:size(y))
       velocity(i,j,:) = [10*sin(x(i))*cos(y(j)), -10*cos(x(i))*sin(y(j))]
     end do
@@ -41,7 +41,7 @@ program advection_diffusion_2D
   procedure(scalar_2D_initializer_i), pointer :: scalar_2D_initializer
   procedure(vector_2D_initializer_i), pointer :: velocity_2D_initializer
   type(scalar_2D_t) s
-  double precision, parameter :: pi = acos(-1D0)
+  real, parameter :: pi = acos(-1E0)
 
   scalar_2D_initializer => scalar_field
   velocity_2D_initializer => taylor_green_velocity
@@ -60,7 +60,7 @@ program advection_diffusion_2D
 
     advance_time: &
     block
-      double precision :: dt = 1D-4
+      real :: dt = 1E-4
       integer step
 
       do step = 1, 500
@@ -88,7 +88,7 @@ pure function d_dt(s, v) result(ds_dt)
     type(scalar_2D_t), intent(in) :: s
     type(vector_2D_t), intent(in) :: v
     type(scalar_2D_t) ds_dt
-    double precision, parameter :: D = 0.5D0
+    real, parameter :: D = 0.5E0
     ds_dt = .div. (D * .grad. s) - .div. (v * s)
   end function
 
diff --git a/example/2D-sink.F90 b/example/2D-sink.F90
index 51a66fa..746279a 100644
--- a/example/2D-sink.F90
+++ b/example/2D-sink.F90
@@ -8,27 +8,27 @@ module sink_2D_functions_m
   implicit none
 
   integer, parameter :: space_dimension = 2
-  double precision, parameter :: pi = acos(-1D0)
+  real, parameter :: pi = acos(-1E0)
 
 contains
 
   pure function velocity(x,y) result(v)
-    double precision, intent(in) :: x(:), y(:)
-    double precision v(size(x),size(y),space_dimension)
-    double precision, parameter :: Q = 1D0
+    real, intent(in) :: x(:), y(:)
+    real v(size(x),size(y),space_dimension)
+    real, parameter :: Q = 1E0
     do concurrent(integer :: i=1:size(x), j=1:size(y))
       associate(r => sqrt(x(i)**2 + y(j)**2))
-        call_julienne_assert(r /= 0D0)
+        call_julienne_assert(r /= 0E0)
         v(i,j,:) = -(Q/(2*pi))*[x(i), y(j)]/(x(i)**2 + y(j)**2)
       end associate
     end do
   end function
 
   pure function divergence(x,y) result(div_v)
-    double precision, intent(in) :: x(:), y(:)
-    double precision div_v(size(x),size(y))
-    call_julienne_assert(.not. any(x == 0D0 .and. y == 0D0))
-    div_v = 0D0
+    real, intent(in) :: x(:), y(:)
+    real div_v(size(x),size(y))
+    call_julienne_assert(.not. any(x == 0E0 .and. y == 0E0))
+    div_v = 0E0
   end function
 
 end module
@@ -45,7 +45,7 @@ program sink_2D
   divergence_2D_initializer => divergence
   vector_2D_initializer => velocity
 
-  associate(v => vector_2D_t(vector_2D_initializer, order=order, cells=[11,11], x_min=[-1D0,-1D0], x_max=[1D0,1D0]))
+  associate(v => vector_2D_t(vector_2D_initializer, order=order, cells=[11,11], x_min=[-1E0,-1E0], x_max=[1E0,1E0]))
     associate(div_v => .div. v, expected_divergence => divergence_2D_t(divergence_2D_initializer, mold=v))
       associate(v_file => v%to_file("v"),div_v_file => div_v%to_file(".div. v"), expected_divergence_file => expected_divergence%to_file("expected .div. v"))
         call v_file%write_lines("example/scripts/sink-velocity.csv") 
diff --git a/example/2D-stagnation-point.F90 b/example/2D-stagnation-point.F90
index d9c5836..5337413 100644
--- a/example/2D-stagnation-point.F90
+++ b/example/2D-stagnation-point.F90
@@ -9,16 +9,16 @@ module velocity_potential_m
 contains
 
   pure function potential(x,y) result(phi)
-    double precision, intent(in) :: x(:), y(:)
-    double precision phi(size(x),size(y))
+    real, intent(in) :: x(:), y(:)
+    real phi(size(x),size(y))
     do concurrent(integer :: j=1:size(y)) default(none) shared(x,y,phi)
       phi(:,j) = (x**2 - y(j)**2)/2
     end do
   end function
 
   pure function potential_gradient(x,y) result(grad_phi)
-    double precision, intent(in) :: x(:), y(:)
-    double precision grad_phi(size(x),size(y),space_dimension)
+    real, intent(in) :: x(:), y(:)
+    real grad_phi(size(x),size(y),space_dimension)
     do concurrent(integer :: i=1:size(x), j=1:size(y))
       grad_phi(i,j,:) = [x(i), -y(j)]
     end do
@@ -38,7 +38,7 @@ program stagnation_point_2D
   vector_2D_initializer => potential_gradient
 
   associate( &
-    phi => scalar_2D_t(scalar_2D_initializer, order=4, cells=[20,20], x_min=[-2D0,-2D0], x_max=[2D0,2D0]) &
+    phi => scalar_2D_t(scalar_2D_initializer, order=4, cells=[20,20], x_min=[-2E0,-2E0], x_max=[2E0,2E0]) &
   )
     associate( &
        v => .grad. phi &
diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90
index 59e8b57..68a501a 100644
--- a/example/burgers-1D.F90
+++ b/example/burgers-1D.F90
@@ -10,8 +10,8 @@ module initial_condition_m
 
   pure function initial_condition(x)
     !! Initial solution to Burgers equation
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: initial_condition(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: initial_condition(:)
     initial_condition = 10*sin(x)
     ! To change this function, please edit only the right-hand-side (RHS) expression,
     ! keeping the rest in place for proper display of the function at runtime.
@@ -66,14 +66,14 @@ program burgers_1D
 
   block
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer
-    double precision, parameter :: pi = acos(-1D0), nu=1D0, t_final=0.6D0
-    double precision, allocatable :: u_surface(:,:)
-    double precision dt
+    real, parameter :: pi = acos(-1E0), nu=1E0, t_final=0.6E0
+    real, allocatable :: u_surface(:,:)
+    real dt
     type(scalar_1D_t) u
     integer step, n
 
     scalar_1D_initializer => initial_condition
-    u = scalar_1D_t(scalar_1D_initializer, order, x_min=0D0, x_max=2*pi, cells=199)
+    u = scalar_1D_t(scalar_1D_initializer, order, x_min=0E0, x_max=2*pi, cells=199)
     dt = diffusion_stability_limit(nu, u%dx(), order)
 
     associate(steps => ceiling(t_final/dt))
@@ -148,17 +148,17 @@ program burgers_1D
 
   pure function d_dt(u, nu) result(du_dt)
     type(scalar_1D_t), intent(in) :: u
-    double precision, intent(in) :: nu
+    real, intent(in) :: nu
     type(scalar_1D_t) du_dt
     du_dt = nu*d2_dx2(u) - d_dx((u**2)/2)
   end function
 
   pure function diffusion_stability_limit(diffusivity,delta_x,order_of_accuracy)  result(stable_time_step)
-    double precision, intent(in) :: diffusivity, delta_x
+    real, intent(in) :: diffusivity, delta_x
     integer, intent(in) :: order_of_accuracy
-    double precision stable_time_step
-    double precision, parameter, dimension(*) :: stability_limit=[2.,2.,2.5,2.79] ! third value needs to be checked
-    double precision, parameter :: safety_factor = 0.9
+    real stable_time_step
+    real, parameter, dimension(*) :: stability_limit=[2.,2.,2.5,2.79] ! third value needs to be checked
+    real, parameter :: safety_factor = 0.9
     ! See Moin, P. (2010) Fundamentals of Engineering Numerical Analysis, 2nd ed., pp. 111-116.
     stable_time_step = safety_factor*stability_limit(order_of_accuracy)*(delta_x**2)/(4*diffusivity)
   end function
diff --git a/example/div-grad-laplacian-1D.F90 b/example/div-grad-laplacian-1D.F90
index 2693631..6d8024a 100644
--- a/example/div-grad-laplacian-1D.F90
+++ b/example/div-grad-laplacian-1D.F90
@@ -8,25 +8,25 @@ module functions_m
 
   pure function f(x)
     !! Function for differentiation
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: f(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: f(:)
     f = (x**3)/6 + (x**2)/2 + 1
     ! To change functions, edit only the right-hand-side (RHS) expression. 
     ! Please keep the rest place for proper display of the function at runtime.
   end function
 
-  double precision elemental function df_dx(x)
+  real elemental function df_dx(x)
     !! 1st-derivative function
-    double precision, intent(in) :: x
+    real, intent(in) :: x
     df_dx = (x**2)/2 + x
     ! To change derivative functions, edit only the RHS.
     ! Please keep the rest place for proper display of the function at runtime.
     ! Also, ensure the new RHS expresses the first derivative of f above.
   end function
 
-  double precision elemental function d2f_dx2(x)
+  real elemental function d2f_dx2(x)
     !! 2nd-derivative function
-    double precision, intent(in) :: x
+    real, intent(in) :: x
     d2f_dx2 = x + 1
     ! To edit the above function, edit only the right-hand side expression,
     ! Please keep the rest place for proper display of the function at runtime.
@@ -93,7 +93,7 @@ program div_grad_laplacian_1D
   subroutine output(order)
     integer, intent(in) :: order
   
-    associate(   s           => scalar_1D_t(scalar_1D_initializer, order=order, cells=20, x_min=0D0, x_max=20D0))
+    associate(   s           => scalar_1D_t(scalar_1D_initializer, order=order, cells=20, x_min=0E0, x_max=20E0))
       associate( grad_s      => .grad. s &
                 ,laplacian_s => .laplacian. s)
         associate( s_grid           => s%grid()      &
@@ -136,7 +136,7 @@ subroutine output(order)
   end subroutine
 
   pure function tabulate(headings, abscissa, expected, actual) result(file)
-    double precision, intent(in), dimension(:) :: abscissa, expected, actual
+    real, intent(in), dimension(:) :: abscissa, expected, actual
     type(string_t), intent(in) :: headings(:)
     type(file_t) file
     integer line
diff --git a/example/extended-gauss-divergence.F90 b/example/extended-gauss-divergence.F90
index 5a57c54..aef7aa1 100644
--- a/example/extended-gauss-divergence.F90
+++ b/example/extended-gauss-divergence.F90
@@ -6,14 +6,14 @@ module integrand_operands_m
 contains
 
   pure function scalar(x) result(f)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: f(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: f(:)
     f = (x**2)/2 ! <-- scalar function
   end function
 
   pure function vector(x) result(v)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: v(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: v(:)
     v = x        ! <-- vector function
   end function
 
@@ -36,7 +36,7 @@ program extended_gauss_divergence
     !! Define default initializations that can be overridden with the command-line arguments
     !! detailed by the usage information below
     integer :: cells_=200, order_=4
-    double precision :: x_min_=0D0, x_max_=1D0
+    real :: x_min_=0E0, x_max_=1E0
   end type
 
   type text_flags_t
@@ -44,7 +44,7 @@ program extended_gauss_divergence
   end type
 
   type(command_line_t) command_line
-  double precision SSS_v_dot_grad_f_dV, SSS_f_div_v_dV, SS_f_v_dot_dA
+  real SSS_v_dot_grad_f_dV, SSS_f_div_v_dV, SS_f_v_dot_dA
 
   if (command_line%argument_present([character(len=len("--help")) :: ("--help"), "-h"])) then
     stop                              new_line('') // new_line('') &
@@ -53,7 +53,7 @@ program extended_gauss_divergence
       // '  --example extended-gauss-divergence \' // new_line('') &
       // '  --compiler flang-new \'                // new_line('') &
       // '  --flag "-O3" \'                        // new_line('') &
-      // '  -- [--help|-h] | [[--cells <integer>] [--order <integer>] [--xmin <double precision>] [--xmax <double precision>] [--div|d] [--grad|g] [--vf|f]]' &
+      // '  -- [--help|-h] | [[--cells <integer>] [--order <integer>] [--xmin <real>] [--xmax <real>] [--div|d] [--grad|g] [--vf|f]]' &
       //                              new_line('') // new_line('') &
       // 'where pipes (|) separate square-bracketed optional arguments and angular brackets indicate user input values.' // new_line('')
   end if
diff --git a/example/print-assembled-1D-operators.F90 b/example/print-assembled-1D-operators.F90
index 9126269..1b217e7 100644
--- a/example/print-assembled-1D-operators.F90
+++ b/example/print-assembled-1D-operators.F90
@@ -33,10 +33,10 @@ program print_assembled_1D_operators
     default_usage: &
     associate(print_all => .not. any([gradient, divergence, len(order)/=0]))
 
-      if (print_all .or. (gradient   .and. len(order)==0) .or. (gradient   .and. order=="2")) call print_gradient_operator(  k=2, dx=1D0, m=16)
-      if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="2")) call print_divergence_operator(k=2, dx=1D0, m=16)
-      if (print_all .or. (gradient   .and. len(order)==0) .or. (gradient   .and. order=="4")) call print_gradient_operator(  k=4, dx=1D0, m=16)
-      if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="4")) call print_divergence_operator(k=4, dx=1D0, m=16)
+      if (print_all .or. (gradient   .and. len(order)==0) .or. (gradient   .and. order=="2")) call print_gradient_operator(  k=2, dx=1E0, m=16)
+      if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="2")) call print_divergence_operator(k=2, dx=1E0, m=16)
+      if (print_all .or. (gradient   .and. len(order)==0) .or. (gradient   .and. order=="4")) call print_gradient_operator(  k=4, dx=1E0, m=16)
+      if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="4")) call print_divergence_operator(k=4, dx=1E0, m=16)
 
     end associate default_usage
   end associate command_line_settings
@@ -45,7 +45,7 @@ program print_assembled_1D_operators
 
   subroutine print_gradient_operator(k, dx, m)
     integer, intent(in) :: k, m 
-    double precision, intent(in) :: dx
+    real, intent(in) :: dx
     integer row
 
     print *, new_line(""), "Gradient operator: order = ", k, " | cells = ", m, " | dx = ", dx
@@ -64,7 +64,7 @@ subroutine print_gradient_operator(k, dx, m)
 
   subroutine print_divergence_operator(k, dx, m)
     integer, intent(in) :: k, m 
-    double precision, intent(in) :: dx
+    real, intent(in) :: dx
     integer row
 
     print *, new_line(""), "Divergence operator: order = ", k, " | cells = ", m, " | dx = ", dx
diff --git a/src/formal/differential_operators_1D_m.F90 b/src/formal/differential_operators_1D_m.F90
index 0b77038..b79f3bf 100644
--- a/src/formal/differential_operators_1D_m.F90
+++ b/src/formal/differential_operators_1D_m.F90
@@ -16,9 +16,9 @@ module differential_operators_1D_m
   type differential_operator_matrix_1D_t
     !! Encapsulate a mimetic matrix
     private
-    double precision, allocatable :: upper_(:,:) !! A  submatrix block (cf. Corbino & Castillo, 2020)
-    double precision, allocatable :: inner_(:)   !! M  submatrix row   (cf. Corbino & Castillo, 2020)
-    double precision, allocatable :: lower_(:,:) !! A' submatrix block (cf. Corbino & Castillo, 2020)
+    real, allocatable :: upper_(:,:) !! A  submatrix block (cf. Corbino & Castillo, 2020)
+    real, allocatable :: inner_(:)   !! M  submatrix row   (cf. Corbino & Castillo, 2020)
+    real, allocatable :: lower_(:,:) !! A' submatrix block (cf. Corbino & Castillo, 2020)
   contains
     procedure, non_overridable :: to_file_t
   end type
@@ -28,9 +28,9 @@ module differential_operators_1D_m
     pure module function construct_matrix_operator(upper, inner, lower) result(differential_operator_matrix_1D)
       !! Construct discrete operator from matrix blocks
       implicit none
-      double precision, intent(in) :: upper(:,:) !! A  submatrix block (cf. Corbino & Castillo, 2020)
-      double precision, intent(in) :: inner(:)   !! M  submatrix row   (cf. Corbino & Castillo, 2020)
-      double precision, intent(in) :: lower(:,:) !! A' submatrix block (cf. Corbino & Castillo, 2020)
+      real, intent(in) :: upper(:,:) !! A  submatrix block (cf. Corbino & Castillo, 2020)
+      real, intent(in) :: inner(:)   !! M  submatrix row   (cf. Corbino & Castillo, 2020)
+      real, intent(in) :: lower(:,:) !! A' submatrix block (cf. Corbino & Castillo, 2020)
       type(differential_operator_matrix_1D_t) differential_operator_matrix_1D
     end function
 
@@ -41,7 +41,7 @@ pure module function construct_matrix_operator(upper, inner, lower) result(diffe
     private
     integer k_ !! order of accuracy
     integer m_ !! number of cells
-    double precision dx_ !! cell width
+    real dx_ !! cell width
   contains
     generic :: operator(.x.) => gradient_matrix_multiply
     procedure, non_overridable, private :: gradient_matrix_multiply
@@ -55,7 +55,7 @@ elemental module function construct_1D_gradient_operator(k, dx, cells) result(gr
       !! Construct a mimetic gradient operator
       implicit none
       integer, intent(in) :: k !! order of accuracy
-      double precision, intent(in) :: dx !! step size
+      real, intent(in) :: dx !! step size
       integer, intent(in) :: cells !! number of grid cells
       type(gradient_operator_1D_t) gradient_operator_1D
     end function
@@ -66,7 +66,7 @@ elemental module function construct_1D_gradient_operator(k, dx, cells) result(gr
     !! Encapsulate kth-order mimetic divergence operator on m_ cells of width dx
     private
     integer k_, m_
-    double precision dx_
+    real dx_
   contains
     generic :: operator(.x.) => divergence_matrix_multiply
     procedure, non_overridable, private :: divergence_matrix_multiply
@@ -81,7 +81,7 @@ elemental module function construct_1D_divergence_operator(k, dx, cells) result(
       !! Construct a mimetic gradient operator
       implicit none
       integer, intent(in) :: k !! order of accuracy
-      double precision, intent(in) :: dx !! step size
+      real, intent(in) :: dx !! step size
       integer, intent(in) :: cells !! number of grid cells
       type(divergence_operator_1D_t) divergence_operator_1D
     end function
@@ -101,30 +101,30 @@ pure module function gradient_matrix_multiply(self, vec) result(matvec_product)
       !! Result is mimetic gradient vector
       implicit none
       class(gradient_operator_1D_t), intent(in) :: self
-      double precision, intent(in) :: vec(:)
-      double precision, allocatable :: matvec_product(:)
+      real, intent(in) :: vec(:)
+      real, allocatable :: matvec_product(:)
     end function
 
     pure module function assemble_gradient(self) result(G)
       !! Result is the assembled 1D mimetic gradient operator matrix
        implicit none
        class(gradient_operator_1D_t), intent(in) :: self
-       double precision, allocatable :: G(:,:)
+       real, allocatable :: G(:,:)
     end function
 
     pure module function assemble_divergence(self) result(D)
       !! Result is the assembled 1D mimetic divergence operator matrix
        implicit none
        class(divergence_operator_1D_t), intent(in) :: self
-       double precision, allocatable :: D(:,:)
+       real, allocatable :: D(:,:)
      end function
 
     pure module function divergence_matrix_multiply(self, vec) result(matvec_product)
       !! Result is mimetic divergence defined at cell centers
       implicit none
       class(divergence_operator_1D_t), intent(in) :: self
-      double precision, intent(in) :: vec(:)
-      double precision, allocatable :: matvec_product(:)
+      real, intent(in) :: vec(:)
+      real, allocatable :: matvec_product(:)
     end function
 
      pure module function to_file_t(self) result(file)
@@ -141,8 +141,8 @@ pure module function to_file_t(self) result(file)
 
   pure function negate_and_flip(A) result(Ap)
     !! Transform a mimetic matrix upper block into a lower block
-    double precision, intent(in) :: A(:,:)
-    double precision, allocatable :: Ap(:,:)
+    real, intent(in) :: A(:,:)
+    real, allocatable :: Ap(:,:)
 
     allocate(Ap, mold=A)
 
diff --git a/src/formal/divergence_1D_s.F90 b/src/formal/divergence_1D_s.F90
index bb52f6f..08313a6 100644
--- a/src/formal/divergence_1D_s.F90
+++ b/src/formal/divergence_1D_s.F90
@@ -40,19 +40,19 @@
   module procedure divergence_1D_weights
       integer c 
 
-      double precision, allocatable :: skin(:)
+      real, allocatable :: skin(:)
 
       select case(self%order_)
       case(2)
-        skin = [double precision::]
+        skin = [real::]
       case(4)
-        skin = [1D0, 2186/1943D0, 1992/2651D0, 1993/1715D0, 649/674D0, 699/700D0, 18170/18171D0, 471744/471745D0]
+        skin = [1E0, 2186/1943E0, 1992/2651E0, 1993/1715E0, 649/674E0, 699/700E0, 18170/18171E0, 471744/471745E0]
       case default
         error stop "unsupported order"
       end select
 
       associate(depth => size(skin))
-        weights = [skin, [(1D0, c = depth+1, self%cells_+2-depth)], skin(depth:1:-1) ] ! m+2 values, where m = self%cells_
+        weights = [skin, [(1E0, c = depth+1, self%cells_+2-depth)], skin(depth:1:-1) ] ! m+2 values, where m = self%cells_
       end associate                                                                    ! cf. Corbino & Castillo (2020) Eqs. 14-15 & 19
 
       call_julienne_assert(self%cells_ .isAtLeast. 2*size(skin))
diff --git a/src/formal/divergence_2D_s.F90 b/src/formal/divergence_2D_s.F90
index ee9abe0..e90aa6c 100644
--- a/src/formal/divergence_2D_s.F90
+++ b/src/formal/divergence_2D_s.F90
@@ -27,7 +27,7 @@
 
     construct_prototype: &
     associate(divergence_1D => divergence_1D_t( &
-       constant = 0D0 &
+       constant = 0E0 &
       ,cells = self%cells_(direction) &
       ,x_min = self%x_min_(direction) &
       ,x_max = self%x_max_(direction) &
diff --git a/src/formal/divergence_3D_s.F90 b/src/formal/divergence_3D_s.F90
index 24d57fc..4aa0d2b 100644
--- a/src/formal/divergence_3D_s.F90
+++ b/src/formal/divergence_3D_s.F90
@@ -28,7 +28,7 @@
 
     construct_prototype: &
     associate(divergence_1D => divergence_1D_t( &
-       constant = 0D0 &
+       constant = 0E0 &
       ,cells = self%cells_(direction) &
       ,x_min = self%x_min_(direction) &
       ,x_max = self%x_max_(direction) &
diff --git a/src/formal/divergence_operator_1D_s.F90 b/src/formal/divergence_operator_1D_s.F90
index 68c23d4..c9cf6dd 100644
--- a/src/formal/divergence_operator_1D_s.F90
+++ b/src/formal/divergence_operator_1D_s.F90
@@ -16,8 +16,8 @@
 
   pure function negate_and_flip(A) result(Ap)
     !! Transform a mimetic matrix upper block into a lower block
-    double precision, intent(in) :: A(:,:)
-    double precision, allocatable :: Ap(:,:)
+    real, intent(in) :: A(:,:)
+    real, allocatable :: Ap(:,:)
     integer row, column
 
     allocate(Ap, mold=A)
@@ -38,7 +38,7 @@ pure function negate_and_flip(A) result(Ap)
  
   module procedure construct_1D_divergence_operator
 
-    double precision, allocatable :: Ap(:,:)
+    real, allocatable :: Ap(:,:)
 
     call_julienne_assert(cells .isAtLeast. 2*k+1)
 
@@ -59,17 +59,17 @@ pure function negate_and_flip(A) result(Ap)
     pure function A_block(k, dx) result(matrix_block)
       !! Compute the upper block submatrix "A" of the Corbino & Castillo (2020) mimetic divergence operator
       integer, intent(in) :: k
-      double precision, intent(in) :: dx
-      double precision, allocatable :: matrix_block(:,:)
+      real, intent(in) :: dx
+      real, allocatable :: matrix_block(:,:)
 
       order_of_accuracy: &
       select case(k)
       case(2)
-        matrix_block = reshape([ double precision :: &
+        matrix_block = reshape([ real :: &
         ], shape=[0,0])
       case(4)
         matrix_block = reshape([ &
-          -11/12D0, 17/24D0, 3/8D0, -5/24D0,  1/24D0 &
+          -11/12E0, 17/24E0, 3/8E0, -5/24E0,  1/24E0 &
         ], shape=[1,5], order=[2,1]) / dx
       case default
         associate(string_k => string_t(k))
@@ -82,15 +82,15 @@ pure function A_block(k, dx) result(matrix_block)
     pure function M(k, dx) result(row)
       !! Compute the middle block submatrix "M" of the Corbino & Castillo (2020) mimetic divergence operator
       integer, intent(in) :: k
-      double precision, intent(in) :: dx
-      double precision, allocatable :: row(:)
+      real, intent(in) :: dx
+      real, allocatable :: row(:)
 
       order_of_accuracy: &
       select case(k)
       case(2)
-        row = [-1D0, 1D0]/ dx        
+        row = [-1E0, 1E0]/ dx        
       case(4)
-        row = [1D0/24D0, -9D0/8D0, 9D0/8D0, -1D0/24D0] / dx        
+        row = [1E0/24E0, -9E0/8E0, 9E0/8E0, -1E0/24E0] / dx        
       case default
         associate(string_k => string_t(k))
           error stop "M (divergence_operator_1D_s): unsupported order of accuracy: " // string_k%string()
@@ -108,7 +108,7 @@ pure function M(k, dx) result(row)
 
   module procedure divergence_matrix_multiply
 
-    double precision, allocatable :: product_inner(:)
+    real, allocatable :: product_inner(:)
 
     associate( &
        upper_rows => size(self%upper_,1) &
@@ -142,11 +142,11 @@ pure function M(k, dx) result(row)
       ,lower_columns => size(self%lower_,2) &
     )
       matvec_product = [ &
-         0D0 &
+         0E0 &
         ,matmul(self%upper_, vec(1 : upper_columns )) &
         ,product_inner &
         ,matmul(self%lower_, vec(size(vec) - lower_columns + 1 : )) &
-        ,0D0 &
+        ,0E0 &
       ]
       call_julienne_assert(size(matvec_product) .equalsExpected. self%m_ + 2)
     end associate
@@ -178,10 +178,10 @@ pure function M(k, dx) result(row)
     pure function e(dir, length) result(unit_vector)
       !! Result is the dir-th column of the len x len identity matrix
       integer, intent(in) :: dir, length
-      double precision :: unit_vector(length)
-      unit_vector(1:dir-1) = 0D0
-      unit_vector(dir)     = 1D0
-      unit_vector(dir+1:)  = 0D0
+      real :: unit_vector(length)
+      unit_vector(1:dir-1) = 0E0
+      unit_vector(dir)     = 1E0
+      unit_vector(dir+1:)  = 0E0
     end function
 
   end procedure
diff --git a/src/formal/gradient_1D_s.F90 b/src/formal/gradient_1D_s.F90
index 010a0ef..b1b018a 100644
--- a/src/formal/gradient_1D_s.F90
+++ b/src/formal/gradient_1D_s.F90
@@ -12,27 +12,27 @@
     ,operator(.within.)
   implicit none
 
-  double precision, parameter :: double_equivalence = 1D-15
+  real, parameter :: real_equivalence = 1E-5
 
 contains
 
   module procedure gradient_1D_weights
 
     integer face
-    double precision, allocatable :: skin(:)
+    real, allocatable :: skin(:)
 
     select case(self%order_)
     case(2)
-      skin = [3/8D0, 9/8D0]
+      skin = [3/8E0, 9/8E0]
     case(4)
-      skin = [227/641D0, 941/766D0, 811/903D0, 1373/1348D0, 1401/1400D0, 36343/36342D0, 943491/943490D0]
+      skin = [227/641E0, 941/766E0, 811/903E0, 1373/1348E0, 1401/1400E0, 36343/36342E0, 943491/943490E0]
     case default
       error stop "unsupported order"
     end select
 
     associate(depth => size(skin))
       call_julienne_assert(self%cells_ .isAtLeast. 2*depth)
-      weights = [skin, [(1D0, face = 1, self%cells_ + 1 - 2*depth)], skin(depth:1:-1) ]
+      weights = [skin, [(1E0, face = 1, self%cells_ + 1 - 2*depth)], skin(depth:1:-1) ]
     end associate
 
     call_julienne_assert(size(weights) .equalsExpected. self%cells_ + 1)
@@ -43,8 +43,8 @@
     call_julienne_assert(size(gradient_1D%values_) .equalsExpected. size(vector_1D%values_))
     call_julienne_assert(gradient_1D%order_ .equalsExpected. vector_1D%order_)
     call_julienne_assert(gradient_1D%cells_ .equalsExpected. vector_1D%cells_)
-    call_julienne_assert(gradient_1D%x_min_ .approximates.    vector_1D%x_min_ .within. double_equivalence)
-    call_julienne_assert(gradient_1D%x_max_ .approximates.    vector_1D%x_max_ .within. double_equivalence)
+    call_julienne_assert(gradient_1D%x_min_ .approximates.    vector_1D%x_min_ .within. real_equivalence)
+    call_julienne_assert(gradient_1D%x_max_ .approximates.    vector_1D%x_max_ .within. real_equivalence)
 
     vector_dot_gradient_1D%tensor_1D_t = tensor_1D_t(   &
        values = gradient_1D%values_ * vector_1D%values_ &
diff --git a/src/formal/gradient_operator_1D_s.F90 b/src/formal/gradient_operator_1D_s.F90
index e758b01..673f50b 100644
--- a/src/formal/gradient_operator_1D_s.F90
+++ b/src/formal/gradient_operator_1D_s.F90
@@ -17,8 +17,8 @@
 
   pure function negate_and_flip(A) result(Ap)
     !! Transform a mimetic matrix upper block into a lower block
-    double precision, intent(in) :: A(:,:)
-    double precision, allocatable :: Ap(:,:)
+    real, intent(in) :: A(:,:)
+    real, allocatable :: Ap(:,:)
     integer row, column
 
     allocate(Ap, mold=A)
@@ -52,17 +52,17 @@ pure function negate_and_flip(A) result(Ap)
 
     pure function corbino_castillo_A(k, dx) result(matrix_block)
       integer, intent(in) :: k
-      double precision, intent(in) :: dx
-      double precision, allocatable :: matrix_block(:,:)
+      real, intent(in) :: dx
+      real, allocatable :: matrix_block(:,:)
 
       order_of_accuracy: &
       select case(k)
       case(2)
-        matrix_block = reshape([-8D0/3D0, 3D0, -1D0/3D0] , shape=[1,3]) / dx
+        matrix_block = reshape([-8E0/3E0, 3E0, -1E0/3E0] , shape=[1,3]) / dx
       case(4)
         matrix_block = reshape([ &
-           -352D0/105D0,  35D0/ 8D0, -35D0/24D0, 21D0/40D0, -5D0/ 56D0 &
-          ,  16D0/105D0, -31D0/24D0,  29D0/24D0, -3D0/40D0,  1D0/168D0 &
+           -352E0/105E0,  35E0/ 8E0, -35E0/24E0, 21E0/40E0, -5E0/ 56E0 &
+          ,  16E0/105E0, -31E0/24E0,  29E0/24E0, -3E0/40E0,  1E0/168E0 &
         ], shape=[2,5], order=[2,1]) / dx
       case default
         associate(string_k => string_t(k))
@@ -74,15 +74,15 @@ pure function corbino_castillo_A(k, dx) result(matrix_block)
 
     pure function corbino_castillo_M(k, dx) result(row)
       integer, intent(in) :: k
-      double precision, intent(in) :: dx
-      double precision, allocatable :: row(:)
+      real, intent(in) :: dx
+      real, allocatable :: row(:)
 
       order_of_accuracy: &
       select case(k)
       case(2)
-        row = [-1D0, 1D0]/ dx        
+        row = [-1E0, 1E0]/ dx        
       case(4)
-        row = [1D0/24D0, -9D0/8D0, 9D0/8D0, -1D0/24D0] / dx        
+        row = [1E0/24E0, -9E0/8E0, 9E0/8E0, -1E0/24E0] / dx        
       case default
         associate(string_k => string_t(k))
           error stop "corbino_castillo_A: unsupported order of accuracy: " // string_k%string()
@@ -95,7 +95,7 @@ pure function corbino_castillo_M(k, dx) result(row)
 
   module procedure gradient_matrix_multiply
 
-    double precision, allocatable :: product_inner(:)
+    real, allocatable :: product_inner(:)
 
     associate( &
        upper_rows => size(self%upper_,1) &
@@ -139,7 +139,7 @@ pure function corbino_castillo_M(k, dx) result(row)
 
     associate(rows => self%m_ + 1, cols => self%m_ + 2)
 
-      allocate(G(rows, cols), source = 0D0)
+      allocate(G(rows, cols), source = 0E0)
 
 #if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT
       do concurrent(integer :: col=1:cols) default(none) shared(G, self, cols)
@@ -160,10 +160,10 @@ pure function corbino_castillo_M(k, dx) result(row)
     pure function e(dir, length) result(unit_vector)
       !! Result is the dir-th column of the len x len identity matrix
       integer, intent(in) :: dir, length
-      double precision :: unit_vector(length)
-      unit_vector(1:dir-1) = 0D0
-      unit_vector(dir)     = 1D0
-      unit_vector(dir+1:)  = 0D0
+      real :: unit_vector(length)
+      unit_vector(1:dir-1) = 0E0
+      unit_vector(dir)     = 1E0
+      unit_vector(dir+1:)  = 0E0
     end function
 
   end procedure
diff --git a/src/formal/interpolator_1D_m.F90 b/src/formal/interpolator_1D_m.F90
index 6b7fec4..3e44c5e 100644
--- a/src/formal/interpolator_1D_m.F90
+++ b/src/formal/interpolator_1D_m.F90
@@ -14,11 +14,11 @@ module interpolator_1D_m
     !! Encapsulate a staggered-grid interpolation matrix with a corresponding matrix-vector product operator
     private
     integer order_, cells_, dx_
-    double precision                 first_
-    double precision, allocatable :: upper_(:,:) 
-    double precision, allocatable :: inner_(:)
-    double precision, allocatable :: lower_(:,:)
-    double precision                 final_
+    real                 first_
+    real, allocatable :: upper_(:,:) 
+    real, allocatable :: inner_(:)
+    real, allocatable :: lower_(:,:)
+    real                 final_
   end type
 
   type, extends(interpolator_1D_t) :: centers_to_faces_1D_t
@@ -37,7 +37,7 @@ elemental module function c2f_constructor(order, cells, dx) result(centers_to_fa
       !! Construct centers-to-faces interpolation operator
       implicit none
       integer, intent(in) :: order, cells
-      double precision, intent(in) :: dx
+      real, intent(in) :: dx
       type(centers_to_faces_1D_t) centers_to_faces_1D
     end function
 
@@ -49,7 +49,7 @@ elemental module function f2c_constructor(order, cells, dx) result(faces_to_cent
       !! Construct centers-to-faces interpolation operator
       implicit none
       integer, intent(in) :: order, cells
-      double precision, intent(in) :: dx
+      real, intent(in) :: dx
       type(faces_to_centers_1D_t) faces_to_centers_1D
     end function
 
@@ -61,16 +61,16 @@ pure module function face_values(self, centers_extended) result(faces)
       !! Interpolate cell-centered values to face-centered values
       implicit none
       class(centers_to_faces_1D_t), intent(in) :: self
-      double precision, intent(in) :: centers_extended(:)
-      double precision, allocatable :: faces(:)
+      real, intent(in) :: centers_extended(:)
+      real, allocatable :: faces(:)
     end function
 
     pure module function center_values_extended(self, faces) result(centers_extended)
       !! Interpolate face-centered values to cell-centered values
       implicit none
       class(faces_to_centers_1D_t), intent(in) :: self
-      double precision, intent(in) :: faces(:)
-      double precision, allocatable :: centers_extended(:)
+      real, intent(in) :: faces(:)
+      real, allocatable :: centers_extended(:)
     end function
 
   end interface
diff --git a/src/formal/interpolator_1D_s.F90 b/src/formal/interpolator_1D_s.F90
index a5ebb2d..15a8bc4 100644
--- a/src/formal/interpolator_1D_s.F90
+++ b/src/formal/interpolator_1D_s.F90
@@ -17,17 +17,17 @@
 
     select case(order)
     case(2)
-      centers_to_faces_1D%first_ = (2D0                                 )/2
-      centers_to_faces_1D%upper_ = (reshape([double precision::], [0,3]))/2
-      centers_to_faces_1D%inner_ = ([1D0,1D0]                           )/2
-      centers_to_faces_1D%lower_ = (reshape([double precision::], [0,3]))/2
-      centers_to_faces_1D%final_ = (2D0                                 )/2
+      centers_to_faces_1D%first_ = (2E0                                 )/2
+      centers_to_faces_1D%upper_ = (reshape([real::], [0,3]))/2
+      centers_to_faces_1D%inner_ = ([1E0,1E0]                           )/2
+      centers_to_faces_1D%lower_ = (reshape([real::], [0,3]))/2
+      centers_to_faces_1D%final_ = (2E0                                 )/2
     case(4)
-      centers_to_faces_1D%first_ = 1D0
-      centers_to_faces_1D%upper_ = reshape([-16,  70, 70, -14,   2], [1,5])/112D0
-      centers_to_faces_1D%inner_ =              [ -7,  63, 63,  -7]        /112D0
-      centers_to_faces_1D%lower_ = reshape([  2, -14, 70,  70, -16], [1,5])/112D0
-      centers_to_faces_1D%final_ = 1D0
+      centers_to_faces_1D%first_ = 1E0
+      centers_to_faces_1D%upper_ = reshape([-16,  70, 70, -14,   2], [1,5])/112E0
+      centers_to_faces_1D%inner_ =              [ -7,  63, 63,  -7]        /112E0
+      centers_to_faces_1D%lower_ = reshape([  2, -14, 70,  70, -16], [1,5])/112E0
+      centers_to_faces_1D%final_ = 1E0
     case default
       error stop "c2f_component_constructor: unsupported order"
     end select
@@ -43,17 +43,17 @@
 
     select case(order)
     case(2)
-      faces_to_centers_1D%first_ = (2D0                               )/2
-      faces_to_centers_1D%upper_ = reshape([double precision::], [0,3])/2
-      faces_to_centers_1D%inner_ =                            [1D0,1D0]/2
-      faces_to_centers_1D%lower_ = reshape([double precision::], [0,3])/2
-      faces_to_centers_1D%final_ = (2D0                               )/2
+      faces_to_centers_1D%first_ = (2E0                               )/2
+      faces_to_centers_1D%upper_ = reshape([real::], [0,3])/2
+      faces_to_centers_1D%inner_ =                            [1E0,1E0]/2
+      faces_to_centers_1D%lower_ = reshape([real::], [0,3])/2
+      faces_to_centers_1D%final_ = (2E0                               )/2
     case(4)
-      faces_to_centers_1D%first_ = 1D0
-      faces_to_centers_1D%upper_ = reshape([35, 140, -70, 28, -5], [1,5])/128D0
-      faces_to_centers_1D%inner_ =         [-8,  72,  72, -8]            /128D0
-      faces_to_centers_1D%lower_ = reshape([-5, 28, -70, 140, 35], [1,5])/128D0
-      faces_to_centers_1D%final_ = 1D0
+      faces_to_centers_1D%first_ = 1E0
+      faces_to_centers_1D%upper_ = reshape([35, 140, -70, 28, -5], [1,5])/128E0
+      faces_to_centers_1D%inner_ =         [-8,  72,  72, -8]            /128E0
+      faces_to_centers_1D%lower_ = reshape([-5, 28, -70, 140, 35], [1,5])/128E0
+      faces_to_centers_1D%final_ = 1E0
     case default
       error stop "f2c_component_constructor: unsupported order"
     end select
diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90
index 1d68f60..1fd5309 100644
--- a/src/formal/scalar_1D_s.F90
+++ b/src/formal/scalar_1D_s.F90
@@ -58,7 +58,7 @@ pure logical function conformable(lhs, rhs)
     type(scalar_1D_t), intent(in) :: lhs, rhs
     call_julienne_assert(size(lhs%values_) .equalsExpected. size(rhs%values_))
     call_julienne_assert(.all.([lhs%cells_,lhs%order_] .equalsExpected. [rhs%cells_,rhs%order_]))
-    call_julienne_assert(.all.([lhs%x_min_,lhs%x_max_] .approximates. [rhs%x_min_,rhs%x_max_] .within. 1D-08))
+    call_julienne_assert(.all.([lhs%x_min_,lhs%x_max_] .approximates. [rhs%x_min_,rhs%x_max_] .within. 1E-5))
     conformable = .true.
   end function
 
@@ -83,7 +83,7 @@ pure logical function conformable(lhs, rhs)
       tensor_1D_t(values =  lhs%values_ + rhs%values_, x_min = rhs%x_min_, x_max = rhs%x_max_, cells = rhs%cells_, order = rhs%order_)
   end procedure
 
-  module procedure premultiply_double
+  module procedure premultiply_real
     lhs_x_rhs%gradient_operator_1D_ = rhs%gradient_operator_1D_
     lhs_x_rhs%tensor_1D_t = &
       tensor_1D_t(values = lhs*rhs%values_, x_min = rhs%x_min_, x_max = rhs%x_max_, cells = rhs%cells_, order = rhs%order_)
@@ -95,7 +95,7 @@ pure logical function conformable(lhs, rhs)
       tensor_1D_t(values = lhs*rhs%values_, x_min = rhs%x_min_, x_max = rhs%x_max_, cells = rhs%cells_, order = rhs%order_)
   end procedure
 
-  module procedure postmultiply_double
+  module procedure postmultiply_real
     lhs_x_rhs%gradient_operator_1D_ = lhs%gradient_operator_1D_
     lhs_x_rhs%tensor_1D_t = &
       tensor_1D_t(values = rhs*lhs%values_, x_min = lhs%x_min_, x_max = lhs%x_max_, cells = lhs%cells_, order = lhs%order_)
@@ -136,8 +136,8 @@ pure logical function conformable(lhs, rhs)
       gradient_1D%tensor_1D_t = tensor_1D_t(self%gradient_operator_1D_ .x. self%values_, self%x_min_, self%x_max_, cells=self%cells_, order=self%order_)
       gradient_1D%divergence_operator_1D_ = divergence_operator_1D_t(self%order_, dx, self%cells_)
       check_corbino_castillo_eq_17: &
-      associate(p => gradient_1D%weights(), b => [-1D0, [(0D0, c = 1, self%cells_)], 1D0])
-        call_julienne_assert((.all. (matmul(transpose(self%gradient_operator_1D_%assemble()), p) .approximates. b/dx .within. 2D-3)))
+      associate(p => gradient_1D%weights(), b => [-1E0, [(0E0, c = 1, self%cells_)], 1E0])
+        call_julienne_assert((.all. (matmul(transpose(self%gradient_operator_1D_%assemble()), p) .approximates. b/dx .within. 2E-3)))
       end associate check_corbino_castillo_eq_17
     end associate
 
diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90
index d23227e..c419e66 100644
--- a/src/formal/scalar_2D_s.F90
+++ b/src/formal/scalar_2D_s.F90
@@ -46,7 +46,7 @@
   module procedure scalar_2D_grid
     call_julienne_assert(self%consistent())
     associate(scalar_1D => scalar_1D_t( &
-       constant = 0D0 &
+       constant = 0E0 &
       ,cells = self%cells_(direction) &
       ,x_min = self%x_min_(direction) &
       ,x_max = self%x_max_(direction) &
@@ -56,7 +56,7 @@
     end associate
   end procedure
 
-  module procedure scalar_2D_postmultiply_double
+  module procedure scalar_2D_postmultiply_real
 
     call_julienne_assert(lhs%consistent())
 
@@ -98,7 +98,7 @@
     call_julienne_assert(lhs_x_rhs%consistent())
   end procedure
 
-  module procedure scalar_2D_premultiply_double
+  module procedure scalar_2D_premultiply_real
     lhs_x_rhs =  rhs * lhs
   end procedure
 
@@ -176,8 +176,8 @@
       gradient_2D%divergence_operator_1D_ = divergence_operator_1D_t(self%order_, dx, self%cells_)
 
      !check_corbino_castillo_eq_17: &
-     !associate(p => gradient_1D%weights(), b => [-1D0, [(0D0, c = 1, self%cells_)], 1D0])
-     !  call_julienne_assert((.all. (matmul(transpose(self%gradient_operator_1D_%assemble()), p) .approximates. b/dx .within. 2D-3)))
+     !associate(p => gradient_1D%weights(), b => [-1E0, [(0E0, c = 1, self%cells_)], 1E0])
+     !  call_julienne_assert((.all. (matmul(transpose(self%gradient_operator_1D_%assemble()), p) .approximates. b/dx .within. 2E-3)))
      !end associate check_corbino_castillo_eq_17
     end associate
 
@@ -199,10 +199,10 @@
        ,y_last => size(rhs%points_(1,1,1,1)%values_,y_dir) - 1 &
      )
        lhs%points_(1,1,1,1)%values_(2:x_last-1, 2:y_last-1) = rhs%points_(1,1,1,1)%values_(2:x_last-1, 2:y_last-1) ! internal points
-       lhs%points_(1,1,1,1)%values_(1         ,  :        ) = 0D0 ! x_min boundary
-       lhs%points_(1,1,1,1)%values_(  x_last  ,  :        ) = 0D0 ! x_max boundary
-       lhs%points_(1,1,1,1)%values_( :        , 1         ) = 0D0 ! y_min boundary
-       lhs%points_(1,1,1,1)%values_( :        ,   y_last  ) = 0D0 ! y_max boundary
+       lhs%points_(1,1,1,1)%values_(1         ,  :        ) = 0E0 ! x_min boundary
+       lhs%points_(1,1,1,1)%values_(  x_last  ,  :        ) = 0E0 ! x_max boundary
+       lhs%points_(1,1,1,1)%values_( :        , 1         ) = 0E0 ! y_min boundary
+       lhs%points_(1,1,1,1)%values_( :        ,   y_last  ) = 0E0 ! y_max boundary
      end associate
 
      lhs%cells_ = rhs%cells_
diff --git a/src/formal/scalar_3D_s.F90 b/src/formal/scalar_3D_s.F90
index d40e5ac..0c964fe 100644
--- a/src/formal/scalar_3D_s.F90
+++ b/src/formal/scalar_3D_s.F90
@@ -46,7 +46,7 @@
   module procedure scalar_3D_grid
     call_julienne_assert(self%consistent())
     associate(scalar_1D => scalar_1D_t( &
-       constant = 0D0 &
+       constant = 0E0 &
       ,cells = self%cells_(direction) &
       ,x_min = self%x_min_(direction) &
       ,x_max = self%x_max_(direction) &
@@ -56,7 +56,7 @@
     end associate
   end procedure
 
-  module procedure scalar_3D_postmultiply_double
+  module procedure scalar_3D_postmultiply_real
 
     call_julienne_assert(lhs%consistent())
 
@@ -98,7 +98,7 @@
     call_julienne_assert(lhs_x_rhs%consistent())
   end procedure
 
-  module procedure scalar_3D_premultiply_double
+  module procedure scalar_3D_premultiply_real
     lhs_x_rhs =  rhs * lhs
   end procedure
 
@@ -186,8 +186,8 @@
       gradient_3D%divergence_operator_1D_ = divergence_operator_1D_t(self%order_, dx, self%cells_)
 
      !check_corbino_castillo_eq_17: &
-     !associate(p => gradient_1D%weights(), b => [-1D0, [(0D0, c = 1, self%cells_)], 1D0])
-     !  call_julienne_assert((.all. (matmul(transpose(self%gradient_operator_1D_%assemble()), p) .approximates. b/dx .within. 3D-3)))
+     !associate(p => gradient_1D%weights(), b => [-1E0, [(0E0, c = 1, self%cells_)], 1E0])
+     !  call_julienne_assert((.all. (matmul(transpose(self%gradient_operator_1D_%assemble()), p) .approximates. b/dx .within. 3E-3)))
      !end associate check_corbino_castillo_eq_17
     end associate
 
@@ -210,12 +210,12 @@
        ,z_last => size(rhs%points_(1,1,1,1)%values_,z_dir) - 1 &
      )
        lhs%points_(1,1,1,1)%values_(2:x_last-1, 2:y_last-1, 2:z_last-1) = rhs%points_(1,1,1,1)%values_(2:x_last-1, 2:y_last-1, 2:z_last-1) ! internal points
-       lhs%points_(1,1,1,1)%values_(1     ,  :   ,  :   ) = 0D0 ! x_min boundary
-       lhs%points_(1,1,1,1)%values_(x_last,  :   ,  :   ) = 0D0 ! x_max boundary
-       lhs%points_(1,1,1,1)%values_(  :   ,  1   ,  :   ) = 0D0 ! y_min boundary
-       lhs%points_(1,1,1,1)%values_(  :   ,y_last,  :   ) = 0D0 ! y_max boundary
-       lhs%points_(1,1,1,1)%values_(  :   ,  :   ,  1   ) = 0D0 ! z_min boundary
-       lhs%points_(1,1,1,1)%values_(  :   ,  :   ,z_last) = 0D0 ! z_max boundary
+       lhs%points_(1,1,1,1)%values_(1     ,  :   ,  :   ) = 0E0 ! x_min boundary
+       lhs%points_(1,1,1,1)%values_(x_last,  :   ,  :   ) = 0E0 ! x_max boundary
+       lhs%points_(1,1,1,1)%values_(  :   ,  1   ,  :   ) = 0E0 ! y_min boundary
+       lhs%points_(1,1,1,1)%values_(  :   ,y_last,  :   ) = 0E0 ! y_max boundary
+       lhs%points_(1,1,1,1)%values_(  :   ,  :   ,  1   ) = 0E0 ! z_min boundary
+       lhs%points_(1,1,1,1)%values_(  :   ,  :   ,z_last) = 0E0 ! z_max boundary
      end associate
 
      lhs%cells_ = rhs%cells_
diff --git a/src/formal/scalar_x_divergence_1D_s.F90 b/src/formal/scalar_x_divergence_1D_s.F90
index 93b9b04..a7de4fa 100644
--- a/src/formal/scalar_x_divergence_1D_s.F90
+++ b/src/formal/scalar_x_divergence_1D_s.F90
@@ -11,7 +11,7 @@
 
   module procedure volume_integrate_scalar_x_divergence_1D
     call_julienne_assert(size(integrand%weights_ ) .equalsExpected. size(integrand%values_)+2)
-    integral  = sum(integrand%weights_ * [0D0, integrand%values_, 0D0])
+    integral  = sum(integrand%weights_ * [0E0, integrand%values_, 0E0])
   end procedure
 
 end submodule scalar_x_divergence_1D_s
diff --git a/src/formal/tensor_2D_s.F90 b/src/formal/tensor_2D_s.F90
index b55d9c6..e8f0594 100644
--- a/src/formal/tensor_2D_s.F90
+++ b/src/formal/tensor_2D_s.F90
@@ -76,8 +76,8 @@
     call_julienne_assert(.all. (shape(self%x_max_) .equalsExpected. shape(tensor_2D%x_max_)))
     call_julienne_assert(.all. (self%cells_ .equalsExpected. tensor_2D%cells_))
     call_julienne_assert(.all. (self%order_ .equalsExpected. tensor_2D%order_))
-    call_julienne_assert(.all. (self%x_min_ .approximates.   tensor_2D%x_min_ .within. 0D0))
-    call_julienne_assert(.all. (self%x_max_ .approximates.   tensor_2D%x_max_ .within. 0D0))
+    call_julienne_assert(.all. (self%x_min_ .approximates.   tensor_2D%x_min_ .within. 0E0))
+    call_julienne_assert(.all. (self%x_max_ .approximates.   tensor_2D%x_max_ .within. 0E0))
     conformable = .true.
   end procedure
 
diff --git a/src/formal/tensor_3D_s.F90 b/src/formal/tensor_3D_s.F90
index 82e1a57..f99f358 100644
--- a/src/formal/tensor_3D_s.F90
+++ b/src/formal/tensor_3D_s.F90
@@ -75,8 +75,8 @@
     call_julienne_assert(.all. (shape(self%x_max_) .equalsExpected. shape(tensor_3D%x_max_)))
     call_julienne_assert(.all. (self%cells_ .equalsExpected. tensor_3D%cells_))
     call_julienne_assert(.all. (self%order_ .equalsExpected. tensor_3D%order_))
-    call_julienne_assert(.all. (self%x_min_ .approximates.   tensor_3D%x_min_ .within. 0D0))
-    call_julienne_assert(.all. (self%x_max_ .approximates.   tensor_3D%x_max_ .within. 0D0))
+    call_julienne_assert(.all. (self%x_min_ .approximates.   tensor_3D%x_min_ .within. 0E0))
+    call_julienne_assert(.all. (self%x_max_ .approximates.   tensor_3D%x_max_ .within. 0E0))
     conformable = .true.
   end procedure
 
diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90
index 471d1f6..7555cc8 100644
--- a/src/formal/tensors_1D_m.F90
+++ b/src/formal/tensors_1D_m.F90
@@ -31,15 +31,15 @@ module tensors_1D_m
     pure function scalar_1D_initializer_i(x) result(f)
       !! Sampling function for initializing a scalar_1D_t object
       implicit none
-      double precision, intent(in) :: x(:)
-      double precision, allocatable :: f(:)
+      real, intent(in) :: x(:)
+      real, allocatable :: f(:)
     end function
 
     pure function vector_1D_initializer_i(x) result(v)
       !! Sampling function for initializing a vector_1D_t object
       implicit none
-      double precision, intent(in) :: x(:)
-      double precision, allocatable :: v(:)
+      real, intent(in) :: x(:)
+      real, allocatable :: v(:)
     end function
 
   end interface
@@ -49,11 +49,11 @@ pure function vector_1D_initializer_i(x) result(v)
     !! Child types define the operations supported by each child, including
     !! gradient (.grad.) for scalars and divergence (.div.) for vectors.
     private
-    double precision x_min_ !! domain lower boundary
-    double precision x_max_ !! domain upper boundary
+    real x_min_ !! domain lower boundary
+    real x_max_ !! domain upper boundary
     integer cells_          !! number of grid cells spanning the domain
     integer order_          !! order of accuracy of mimetic discretization
-    double precision, allocatable :: values_(:) !! tensor components at spatial locations
+    real, allocatable :: values_(:) !! tensor components at spatial locations
   contains
     procedure, non_overridable, private :: is_cell_centered
     procedure, non_overridable, private :: is_face_centered
@@ -70,9 +70,9 @@ pure function vector_1D_initializer_i(x) result(v)
     pure module function construct_1D_tensor_from_components(values, x_min, x_max, cells, order) result(tensor_1D)
       !! User-defined constructor: result is a 1D tensor defined by assigning the dummy arguments to corresponding components
       implicit none
-      double precision, intent(in) :: values(:) !! tensor components at grid locations define by child
-      double precision, intent(in) :: x_min     !! grid location minimum
-      double precision, intent(in) :: x_max     !! grid location maximum
+      real, intent(in) :: values(:) !! tensor components at grid locations define by child
+      real, intent(in) :: x_min     !! grid location minimum
+      real, intent(in) :: x_max     !! grid location maximum
       integer,          intent(in) :: cells     !! number of grid cells spanning the domain
       integer,          intent(in) :: order     !! order of accuracy
       type(tensor_1D_t) tensor_1D
@@ -90,7 +90,7 @@ pure module function construct_1D_tensor_from_components(values, x_min, x_max, c
     generic :: operator(-) => subtract_scalar_1D
     generic :: operator(+) => add_scalar_1D
     generic :: operator(/) => divide_by_integer
-    generic :: operator(*) => premultiply_double, postmultiply_double, premultiply_integer, postmultiply_integer, multiply_1D_scalars
+    generic :: operator(*) => premultiply_real, postmultiply_real, premultiply_integer, postmultiply_integer, multiply_1D_scalars
     generic :: operator(**) => exponentiate
     generic :: operator(.grad.) => grad
     generic :: operator(.laplacian.) => laplacian
@@ -100,8 +100,8 @@ pure module function construct_1D_tensor_from_components(values, x_min, x_max, c
     procedure, non_overridable, private :: scalar_1D_grid
     procedure, non_overridable, private :: divide_by_integer 
     procedure, non_overridable, private :: multiply_1D_scalars
-    procedure, non_overridable, private, pass(rhs) :: premultiply_double
-    procedure, non_overridable, private :: postmultiply_double
+    procedure, non_overridable, private, pass(rhs) :: premultiply_real
+    procedure, non_overridable, private :: postmultiply_real
     procedure, non_overridable, private, pass(rhs) :: premultiply_integer
     procedure, non_overridable, private :: postmultiply_integer
     procedure, non_overridable, private :: exponentiate
@@ -117,19 +117,19 @@ pure module function construct_1D_scalar_from_function(initializer, order, cells
       procedure(scalar_1D_initializer_i), pointer :: initializer
       integer, intent(in) :: order !! order of accuracy
       integer, intent(in) :: cells !! number of grid cells spanning the domain
-      double precision, intent(in) :: x_min !! grid location minimum
-      double precision, intent(in) :: x_max !! grid location maximum
+      real, intent(in) :: x_min !! grid location minimum
+      real, intent(in) :: x_max !! grid location maximum
       type(scalar_1D_t) scalar_1D
     end function
 
     pure module function construct_1D_scalar_constant(constant, order, cells, x_min, x_max) result(scalar_1D)
       !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator
       implicit none
-      double precision, intent(in) :: constant !! scalar value
+      real, intent(in) :: constant !! scalar value
       integer, intent(in) :: order !! order of accuracy
       integer, intent(in) :: cells !! number of grid cells spanning the domain
-      double precision, intent(in) :: x_min !! grid location minimum
-      double precision, intent(in) :: x_max !! grid location maximum
+      real, intent(in) :: x_min !! grid location minimum
+      real, intent(in) :: x_max !! grid location maximum
       type(scalar_1D_t) scalar_1D
     end function
 
@@ -177,8 +177,8 @@ pure module function construct_1D_vector_from_function(initializer, order, cells
       procedure(vector_1D_initializer_i), pointer :: initializer
       integer, intent(in) :: order !! order of accuracy
       integer, intent(in) :: cells !! number of grid cells spanning the domain
-      double precision, intent(in) :: x_min !! grid location minimum
-      double precision, intent(in) :: x_max !! grid location maximum
+      real, intent(in) :: x_min !! grid location minimum
+      real, intent(in) :: x_max !! grid location maximum
       type(vector_1D_t) vector_1D
     end function
 
@@ -191,11 +191,11 @@ pure module function construct_1D_vector_from_parent(tensor_1D) result(vector_1D
     pure module function construct_1D_vector_constant(constant, order, cells, x_min, x_max) result(vector_1D)
       !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator
       implicit none
-      double precision, intent(in) :: constant !! scalar value
+      real, intent(in) :: constant !! scalar value
       integer, intent(in) :: order !! order of accuracy
       integer, intent(in) :: cells !! number of grid cells spanning the domain
-      double precision, intent(in) :: x_min !! grid location minimum
-      double precision, intent(in) :: x_max !! grid location maximum
+      real, intent(in) :: x_min !! grid location minimum
+      real, intent(in) :: x_max !! grid location maximum
       type(vector_1D_t) vector_1D
     end function
 
@@ -214,7 +214,7 @@ pure module function construct_1D_vector_constant(constant, order, cells, x_min,
   type, extends(tensor_1D_t) :: vector_dot_gradient_1D_t
     !! Result is the dot product of a 1D vector field and a 1D gradient field
     private
-    double precision, allocatable :: weights_(:)
+    real, allocatable :: weights_(:)
   contains
     generic :: operator(.SSS.) => volume_integrate_vector_dot_grad_scalar_1D
     procedure, non_overridable, private, pass(integrand) ::volume_integrate_vector_dot_grad_scalar_1D
@@ -237,11 +237,11 @@ pure module function construct_1D_vector_constant(constant, order, cells, x_min,
 
     pure module function construct_1D_divergence_constant(constant, order, cells, x_min, x_max) result(divergence_1D)
       implicit none
-      double precision, intent(in) :: constant !! scalar value
+      real, intent(in) :: constant !! scalar value
       integer, intent(in) :: order !! order of accuracy
       integer, intent(in) :: cells !! number of grid cells spanning the domain
-      double precision, intent(in) :: x_min !! grid location minimum
-      double precision, intent(in) :: x_max !! grid location maximum
+      real, intent(in) :: x_min !! grid location minimum
+      real, intent(in) :: x_max !! grid location maximum
       type(divergence_1D_t) divergence_1D
     end function
 
@@ -250,7 +250,7 @@ pure module function construct_1D_divergence_constant(constant, order, cells, x_
   type, extends(tensor_1D_t) :: scalar_x_divergence_1D_t
     !! product of a 1D scalar field and a 1D divergence field
     private
-    double precision, allocatable :: weights_(:)
+    real, allocatable :: weights_(:)
   contains
     generic :: operator(.SSS.) => volume_integrate_scalar_x_divergence_1D
     procedure, non_overridable, private, pass(integrand) :: volume_integrate_scalar_x_divergence_1D
@@ -288,7 +288,7 @@ pure module function dA(self)
       !! .SS. (f .x. (v .dot. dA))
       implicit none
       class(vector_1D_t), intent(in) :: self
-      double precision dA
+      real dA
     end function
 
     pure module function cells(self)
@@ -302,49 +302,49 @@ pure module function dx(self)
       !! Result is the uniform cell width
       implicit none
       class(tensor_1D_t), intent(in) :: self
-      double precision dx
+      real dx
     end function
 
     pure module function scalar_1D_grid(self) result(cell_centers_extended)
       !! Result is the array of locations at which 1D scalars are defined: cell centers agumented by spatial boundaries
       implicit none
       class(scalar_1D_t), intent(in) :: self
-      double precision, allocatable :: cell_centers_extended(:)
+      real, allocatable :: cell_centers_extended(:)
     end function
 
     pure module function vector_1D_grid(self) result(cell_faces)
       !! Result is the array of cell face locations (of unit area for 1D) at which 1D vectors are defined
       implicit none
       class(vector_1D_t), intent(in) :: self
-      double precision, allocatable :: cell_faces(:)
+      real, allocatable :: cell_faces(:)
     end function
 
     pure module function divergence_1D_grid(self) result(cell_centers)
       !! Result is the array of cell centers at which 1D divergences are defined
       implicit none
       class(divergence_1D_t), intent(in) :: self
-      double precision, allocatable :: cell_centers(:)
+      real, allocatable :: cell_centers(:)
     end function
 
     pure module function scalar_1D_values(self) result(cell_centers_extended_values)
       !! Result is an array of 1D scalar values at boundaries and cell centers
       implicit none
       class(scalar_1D_t), intent(in) :: self
-      double precision, allocatable :: cell_centers_extended_values(:)
+      real, allocatable :: cell_centers_extended_values(:)
     end function
 
     pure module function vector_1D_values(self) result(face_centered_values)
       !! Result is an array of the 1D vector values at cell faces (of unit area 1D)
       implicit none
       class(vector_1D_t), intent(in) :: self
-      double precision, allocatable :: face_centered_values(:)
+      real, allocatable :: face_centered_values(:)
     end function
 
     pure module function divergence_1D_values(self) result(cell_centered_values)
       !! Result is an array of 1D divergences at cell centers
       implicit none
       class(divergence_1D_t), intent(in) :: self
-      double precision, allocatable :: cell_centered_values(:)
+      real, allocatable :: cell_centered_values(:)
     end function
 
     pure module function d_dx(self) result(dself_dx)
@@ -368,10 +368,10 @@ pure module function grad(self) result(gradient_1D)
       type(gradient_1D_t) gradient_1D
     end function
 
-    pure module function premultiply_double(lhs, rhs) result(lhs_x_rhs)
+    pure module function premultiply_real(lhs, rhs) result(lhs_x_rhs)
       !! Result is the product of lhs and rhs
       implicit none
-      double precision, intent(in) :: lhs 
+      real, intent(in) :: lhs 
       class(scalar_1D_t), intent(in) :: rhs
       type(scalar_1D_t) lhs_x_rhs
     end function
@@ -384,11 +384,11 @@ pure module function premultiply_integer(lhs, rhs) result(lhs_x_rhs)
       type(scalar_1D_t) lhs_x_rhs
     end function
 
-    pure module function postmultiply_double(lhs, rhs) result(lhs_x_rhs)
+    pure module function postmultiply_real(lhs, rhs) result(lhs_x_rhs)
       !! Result is the product of lhs and rhs
       implicit none
       class(scalar_1D_t), intent(in) :: lhs
-      double precision, intent(in) :: rhs 
+      real, intent(in) :: rhs 
       type(scalar_1D_t) lhs_x_rhs
     end function
 
@@ -465,21 +465,21 @@ pure module function volume_integrate_vector_dot_grad_scalar_1D(integrand) resul
       !! Result is the mimetic quadrature corresponding to a volume integral of a vector-gradient dot product
       implicit none
       class(vector_dot_gradient_1D_t), intent(in) :: integrand
-      double precision integral
+      real integral
     end function
 
     pure module function volume_integrate_scalar_x_divergence_1D(integrand) result(integral)
       !! Result is the mimetic quadrature corresponding to a volume integral of a scalar-divergence product
       implicit none
       class(scalar_x_divergence_1D_t), intent(in) :: integrand
-      double precision integral
+      real integral
     end function
 
     pure module function surface_integrate_vector_x_scalar_1D(integrand) result(integral)
       !! Result is the mimetic quadrature correspondingto a surface integral of a scalar-vector product
       implicit none
       class(weighted_product_1D_t), intent(in) :: integrand
-      double precision integral
+      real integral
     end function
 
     pure module function dot(vector_1D, gradient_1D) result(vector_dot_gradient_1D)
@@ -496,7 +496,7 @@ pure module function dot_surface_normal(vector_1D, dS) result(v_dot_dS)
       !! The sign of the dot-product is incorporated into the weights in the weighted multiplication operator(.x.).
       implicit none
       class(vector_1D_t), intent(in) :: vector_1D
-      double precision, intent(in) :: dS
+      real, intent(in) :: dS
       type(vector_1D_t) v_dot_dS
     end function
 
@@ -513,7 +513,7 @@ pure module function gradient_1D_weights(self) result(weights)
       !! inner product  of a vector_1D_t object and a gradient_1D_t object.
       implicit none
       class(tensor_1D_t), intent(in) :: self
-      double precision, allocatable :: weights(:)
+      real, allocatable :: weights(:)
     end function
 
     pure module function divergence_1D_weights(self) result(weights)
@@ -521,7 +521,7 @@ pure module function divergence_1D_weights(self) result(weights)
       !! inner product  of a vector_1D_t object and a gradient_1D_t object.
       implicit none
       class(tensor_1D_t), intent(in) :: self
-      double precision, allocatable :: weights(:)
+      real, allocatable :: weights(:)
     end function
 
     pure module function premultiply_scalar_1D(scalar_1D, divergence_1D) result(scalar_x_divergence_1D)
@@ -545,16 +545,16 @@ pure module function postmultiply_scalar_1D(divergence_1D, scalar_1D) result(sca
 contains
 
   pure function cell_centers_extended_1D(x_min, x_max, cells) result(x)
-    double precision, intent(in) :: x_min, x_max
+    real, intent(in) :: x_min, x_max
     integer, intent(in) :: cells
-    double precision, allocatable:: x(:)
+    real, allocatable:: x(:)
     x = [x_min, cell_centers_1D(x_min, x_max, cells), x_max]
   end function
 
   pure function faces_1D(x_min, x_max, cells) result(x)
-    double precision, intent(in) :: x_min, x_max
+    real, intent(in) :: x_min, x_max
     integer, intent(in) :: cells
-    double precision, allocatable:: x(:)
+    real, allocatable:: x(:)
     integer cell
     associate(dx => (x_max - x_min)/cells)
       x = [x_min, x_min + [(cell*dx, cell = 1, cells-1)], x_max]
@@ -562,9 +562,9 @@ pure function faces_1D(x_min, x_max, cells) result(x)
   end function
 
   pure function cell_centers_1D(x_min, x_max, cells) result(x)
-    double precision, intent(in) :: x_min, x_max
+    real, intent(in) :: x_min, x_max
     integer, intent(in) :: cells
-    double precision, allocatable:: x(:)
+    real, allocatable:: x(:)
     integer cell
     associate(dx => (x_max - x_min)/cells)
       x = x_min + dx/2. + [((cell-1)*dx, cell = 1, cells)]
diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90
index 8dfeef2..0751710 100644
--- a/src/formal/tensors_2D_m.F90
+++ b/src/formal/tensors_2D_m.F90
@@ -27,29 +27,29 @@ module tensors_2D_m
     pure function scalar_2D_initializer_i(x,y) result(f)
       !! Sampling function for initializing a scalar_2D_t object
       implicit none
-      double precision, intent(in) :: x(:), y(:)
-      double precision f(size(x),size(y))
+      real, intent(in) :: x(:), y(:)
+      real f(size(x),size(y))
     end function
 
     pure function divergence_2D_initializer_i(x,y) result(f)
       !! Sampling function for initializing a divergence_2D_t object
       implicit none
-      double precision, intent(in) :: x(:), y(:)
-      double precision f(size(x),size(y))
+      real, intent(in) :: x(:), y(:)
+      real f(size(x),size(y))
     end function
 
     pure function vector_2D_initializer_i(x,y) result(v)
       !! Sampling function for initializing a vector_2D_t object
       import space_dimension
       implicit none
-      double precision, intent(in) :: x(:), y(:)
-      double precision v(size(x),size(y),space_dimension)
+      real, intent(in) :: x(:), y(:)
+      real v(size(x),size(y),space_dimension)
     end function
 
   end interface
 
   type points_2D_t
-    double precision, allocatable :: values_(:,:) !! tensor component values at 2D locations
+    real, allocatable :: values_(:,:) !! tensor component values at 2D locations
   end type
 
   type tensor_2D_t
@@ -58,8 +58,8 @@ pure function vector_2D_initializer_i(x,y) result(v)
     !! gradient (.grad.) for scalars and divergence (.div.) for vectors.
     private
     type(points_2D_t), allocatable :: points_(:,:,:,:) !! tensor values indexable up to rank 4
-    double precision x_min_(space_dimension) !! domain lower boundary
-    double precision x_max_(space_dimension) !! domain upper boundary
+    real x_min_(space_dimension) !! domain lower boundary
+    real x_max_(space_dimension) !! domain upper boundary
     integer cells_(space_dimension) !! number of grid cells spanning the domain
     integer order_ !! order of accuracy of mimetic discretization
   contains
@@ -74,8 +74,8 @@ pure function vector_2D_initializer_i(x,y) result(v)
     pure module function construct_2D_tensor_from_components(points, cells, x_min, x_max, order) result(tensor_2D)
       implicit none
       type(points_2D_t), intent(in) :: points(:,:,:,:) !! tensor values at 2D spatial locations
-      double precision, intent(in) :: x_min(:) !! domain lower boundary
-      double precision, intent(in) :: x_max(:) !! domain upper boundary
+      real, intent(in) :: x_min(:) !! domain lower boundary
+      real, intent(in) :: x_max(:) !! domain upper boundary
       integer, intent(in) :: cells(:) !! number of grid cells spanning the domain
       integer, intent(in) :: order !! order of accuracy of mimetic discretization
       type(tensor_2D_t) tensor_2D
@@ -90,7 +90,7 @@ pure module function construct_2D_tensor_from_components(points, cells, x_min, x
   contains
     generic :: assignment(=) => scalar_2D_assign_divergence
     generic :: operator(.grad.) => scalar_2D_gradient
-    generic :: operator(*) => scalar_2D_postmultiply_double, scalar_2D_premultiply_double &
+    generic :: operator(*) => scalar_2D_postmultiply_real, scalar_2D_premultiply_real &
                              ,scalar_2D_postmultiply_integer, scalar_2D_premultiply_integer
     generic :: operator(+) => scalar_2D_plus_scalar
     generic :: values => scalar_2D_values
@@ -105,9 +105,9 @@ pure module function construct_2D_tensor_from_components(points, cells, x_min, x
     procedure, non_overridable, private :: scalar_2D_values
     procedure, non_overridable, private :: scalar_2D_grid
     procedure, non_overridable, private :: scalar_2D_consistent
-    procedure, non_overridable, private :: scalar_2D_postmultiply_double, scalar_2D_postmultiply_integer
+    procedure, non_overridable, private :: scalar_2D_postmultiply_real, scalar_2D_postmultiply_integer
     procedure, non_overridable, private :: scalar_2D_plus_scalar
-    procedure, non_overridable, private, pass(rhs) :: scalar_2D_premultiply_double, scalar_2D_premultiply_integer
+    procedure, non_overridable, private, pass(rhs) :: scalar_2D_premultiply_real, scalar_2D_premultiply_integer
   end type
 
   interface scalar_2D_t
@@ -118,8 +118,8 @@ pure module function construct_2D_scalar_from_function(initializer, order, cells
       procedure(scalar_2D_initializer_i), pointer :: initializer
       integer, intent(in) :: order !! order of accuracy
       integer, intent(in) :: cells(:) !! number of grid cells spanning each spatial direction
-      double precision, intent(in) :: x_min(:) !! grid location minima
-      double precision, intent(in) :: x_max(:) !! grid location maxima
+      real, intent(in) :: x_min(:) !! grid location minima
+      real, intent(in) :: x_max(:) !! grid location maxima
       type(scalar_2D_t) scalar_2D
     end function
 
@@ -182,8 +182,8 @@ pure module function construct_2D_vector_from_function(initializer, order, cells
       procedure(vector_2D_initializer_i), pointer :: initializer
       integer, intent(in) :: order !! order of accuracy
       integer, intent(in) :: cells(:) !! number of grid cells spanning each spatial direction
-      double precision, intent(in) :: x_min(:) !! grid location minima
-      double precision, intent(in) :: x_max(:) !! grid location maxima
+      real, intent(in) :: x_min(:) !! grid location minima
+      real, intent(in) :: x_max(:) !! grid location maxima
       type(vector_2D_t) vector_2D
     end function
 
@@ -253,8 +253,8 @@ pure module function construct_2D_divergence_from_function(initializer, order, c
       procedure(scalar_2D_initializer_i), pointer, intent(in) :: initializer
       integer, intent(in) :: order !! order of accuracy
       integer, intent(in) :: cells(:) !! number of grid cells spanning each spatial direction
-      double precision, intent(in) :: x_min(:) !! grid location minima
-      double precision, intent(in) :: x_max(:) !! grid location maxima
+      real, intent(in) :: x_min(:) !! grid location minima
+      real, intent(in) :: x_max(:) !! grid location maxima
       type(divergence_2D_t) divergence_2D
     end function
 
@@ -303,14 +303,14 @@ pure module function scalar_2D_to_faces(self, direction) result(scalars)
       implicit none
       class(scalar_2D_t), intent(in) :: self
       integer, intent(in) :: direction
-      double precision, allocatable :: scalars(:,:)
+      real, allocatable :: scalars(:,:)
     end function
 
     pure module function scalar_2D_values(self) result(values)
       !! Scalar values getter
       implicit none
       class(scalar_2D_t), intent(in) :: self
-      double precision, allocatable :: values(:,:)
+      real, allocatable :: values(:,:)
     end function
 
     pure module function scalar_2D_grid(self, direction) result(scalar_grid_1D)
@@ -318,35 +318,35 @@ pure module function scalar_2D_grid(self, direction) result(scalar_grid_1D)
       implicit none
       class(scalar_2D_t), intent(in) :: self
       integer, intent(in) :: direction
-      double precision, allocatable :: scalar_grid_1D(:)
+      real, allocatable :: scalar_grid_1D(:)
     end function
 
-    pure module function scalar_2D_postmultiply_double(lhs, rhs) result(lhs_x_rhs)
-      !! Result is product of the scalar_2D_t lhs and double-precision (constant) rhs
+    pure module function scalar_2D_postmultiply_real(lhs, rhs) result(lhs_x_rhs)
+      !! Result is product of the scalar_2D_t lhs and real (constant) rhs
       implicit none
       class(scalar_2D_t), intent(in) :: lhs
-      double precision, intent(in) :: rhs
+      real, intent(in) :: rhs
       type(scalar_2D_t) lhs_x_rhs
     end function
 
     pure module function scalar_2D_postmultiply_integer(lhs, rhs) result(lhs_x_rhs)
-      !! Result is product of the scalar_2D_t lhs and double-precision (constant) rhs
+      !! Result is product of the scalar_2D_t lhs and real (constant) rhs
       implicit none
       class(scalar_2D_t), intent(in) :: lhs
       integer, intent(in) :: rhs
       type(scalar_2D_t) lhs_x_rhs
     end function
 
-    pure module function scalar_2D_premultiply_double(lhs, rhs) result(lhs_x_rhs)
-      !! Result is product of the scalar_2D_t rhs and double-precision (constant) lhs
+    pure module function scalar_2D_premultiply_real(lhs, rhs) result(lhs_x_rhs)
+      !! Result is product of the scalar_2D_t rhs and real (constant) lhs
       implicit none
       class(scalar_2D_t), intent(in) :: rhs
-      double precision, intent(in) :: lhs
+      real, intent(in) :: lhs
       type(scalar_2D_t) lhs_x_rhs
     end function
 
     pure module function scalar_2D_premultiply_integer(lhs, rhs) result(lhs_x_rhs)
-      !! Result is product of the scalar_2D_t rhs and double-precision (constant) lhs
+      !! Result is product of the scalar_2D_t rhs and real (constant) lhs
       implicit none
       class(scalar_2D_t), intent(in) :: rhs
       integer, intent(in) :: lhs
@@ -354,7 +354,7 @@ pure module function scalar_2D_premultiply_integer(lhs, rhs) result(lhs_x_rhs)
     end function
 
     pure module function scalar_2D_plus_scalar(lhs, rhs) result(lhs_plus_rhs)
-      !! Result is product of the scalar_2D_t lhs and double-precision (constant) rhs
+      !! Result is product of the scalar_2D_t lhs and real (constant) rhs
       implicit none
       class(scalar_2D_t), intent(in) :: lhs, rhs
       type(scalar_2D_t) lhs_plus_rhs
@@ -365,7 +365,7 @@ pure module function vector_2D_grid(self, component, coordinate) result(vector_g
       implicit none
       class(vector_2D_t), intent(in) :: self
       integer, intent(in) :: component, coordinate
-      double precision, allocatable :: vector_grid_1D(:) !! grid points along the requested coordinate direction
+      real, allocatable :: vector_grid_1D(:) !! grid points along the requested coordinate direction
     end function
 
     pure module function vector_2D_values(self, direction) result(vector_values)
@@ -373,7 +373,7 @@ pure module function vector_2D_values(self, direction) result(vector_values)
       implicit none
       class(vector_2D_t), intent(in) :: self
       integer, intent(in) :: direction
-      double precision, allocatable :: vector_values(:,:) 
+      real, allocatable :: vector_values(:,:) 
     end function
 
     pure module function divergence_2D_grid(self, direction) result(divergence_grid_1D)
@@ -381,14 +381,14 @@ pure module function divergence_2D_grid(self, direction) result(divergence_grid_
       implicit none
       class(divergence_2D_t), intent(in) :: self
       integer, intent(in) :: direction
-      double precision, allocatable :: divergence_grid_1D(:) !! grid points along the requested coordinate direction
+      real, allocatable :: divergence_grid_1D(:) !! grid points along the requested coordinate direction
     end function
 
     pure module function vector_2D_to_centers_extended(self) result(vectors)
       !! Vector values getter
       implicit none
       class(vector_2D_t), intent(in) :: self
-      double precision, allocatable :: vectors(:,:,:)
+      real, allocatable :: vectors(:,:,:)
     end function
 
     pure module function vector_2D_consistent(self) result(self_consistent)
@@ -402,7 +402,7 @@ pure module function divergence_2D_values(self) result(divergences)
       !! Vector values getter
       implicit none
       class(divergence_2D_t), intent(in) :: self
-      double precision, allocatable :: divergences(:,:)
+      real, allocatable :: divergences(:,:)
     end function
 
     pure module function scalar_2D_gradient(self) result(gradient_2D)
@@ -446,7 +446,7 @@ pure module function gradient_2D_postmultiply_constant(lhs, rhs) result(product)
       !! Result is product of the gradient_2D_t lhs and the constant rhs
       implicit none
       class(gradient_2D_t), intent(in) :: lhs
-      double precision, intent(in) :: rhs
+      real, intent(in) :: rhs
       type(gradient_2D_t) product
     end function
 
@@ -454,7 +454,7 @@ pure module function gradient_2D_premultiply_constant(lhs, rhs) result(product)
       !! Result is product of the gradient_2D_t rhs and the constant lhs
       implicit none
       class(gradient_2D_t), intent(in) :: rhs
-      double precision, intent(in) :: lhs
+      real, intent(in) :: lhs
       type(gradient_2D_t) product
     end function
 
@@ -462,7 +462,7 @@ pure module function divergence_2D_postmultiply_constant(lhs, rhs) result(lhs_x_
       !! Result is product of the divergence_2D_t lhs and the constant rhs
       implicit none
       class(divergence_2D_t), intent(in) :: lhs
-      double precision, intent(in) :: rhs
+      real, intent(in) :: rhs
       type(divergence_2D_t) lhs_x_rhs
     end function
 
@@ -470,7 +470,7 @@ pure module function divergence_2D_premultiply_constant(lhs, rhs) result(lhs_x_r
       !! Result is product of the constant rhs and the lhs divergence_2D_t
       implicit none
       class(divergence_2D_t), intent(in) :: rhs
-      double precision, intent(in) :: lhs
+      real, intent(in) :: lhs
       type(divergence_2D_t) lhs_x_rhs
     end function
 
diff --git a/src/formal/tensors_3D_m.F90 b/src/formal/tensors_3D_m.F90
index 41af54f..cdf629f 100644
--- a/src/formal/tensors_3D_m.F90
+++ b/src/formal/tensors_3D_m.F90
@@ -26,29 +26,29 @@ module tensors_3D_m
     pure function scalar_3D_initializer_i(x,y,z) result(f)
       !! Sampling function for initializing a scalar_3D_t object
       implicit none
-      double precision, intent(in) :: x(:), y(:), z(:)
-      double precision f(size(x),size(y),size(z))
+      real, intent(in) :: x(:), y(:), z(:)
+      real f(size(x),size(y),size(z))
     end function
 
     pure function divergence_3D_initializer_i(x,y,z) result(f)
       !! Sampling function for initializing a divergence_3D_t object
       implicit none
-      double precision, intent(in) :: x(:), y(:), z(:)
-      double precision f(size(x),size(y),size(z))
+      real, intent(in) :: x(:), y(:), z(:)
+      real f(size(x),size(y),size(z))
     end function
 
     pure function vector_3D_initializer_i(x,y,z) result(v)
       !! Sampling function for initializing a vector_3D_t object
       import space_dimension
       implicit none
-      double precision, intent(in) :: x(:), y(:), z(:)
-      double precision v(size(x),size(y),size(z),space_dimension)
+      real, intent(in) :: x(:), y(:), z(:)
+      real v(size(x),size(y),size(z),space_dimension)
     end function
 
   end interface
 
   type points_3D_t
-    double precision, allocatable :: values_(:,:,:) !! tensor component values at 3D locations
+    real, allocatable :: values_(:,:,:) !! tensor component values at 3D locations
   end type
 
   type tensor_3D_t
@@ -57,8 +57,8 @@ pure function vector_3D_initializer_i(x,y,z) result(v)
     !! gradient (.grad.) for scalars and divergence (.div.) for vectors.
     private
     type(points_3D_t), allocatable :: points_(:,:,:,:) !! tensor values indexable up to rank 4
-    double precision x_min_(space_dimension) !! domain lower boundary
-    double precision x_max_(space_dimension) !! domain upper boundary
+    real x_min_(space_dimension) !! domain lower boundary
+    real x_max_(space_dimension) !! domain upper boundary
     integer cells_(space_dimension) !! number of grid cells spanning the domain
     integer order_ !! order of accuracy of mimetic discretization
   contains
@@ -73,8 +73,8 @@ pure function vector_3D_initializer_i(x,y,z) result(v)
     pure module function construct_3D_tensor_from_components(points, cells, x_min, x_max, order) result(tensor_3D)
       implicit none
       type(points_3D_t), intent(in) :: points(:,:,:,:) !! tensor values at 3D spatial locations
-      double precision, intent(in) :: x_min(:) !! domain lower boundary
-      double precision, intent(in) :: x_max(:) !! domain upper boundary
+      real, intent(in) :: x_min(:) !! domain lower boundary
+      real, intent(in) :: x_max(:) !! domain upper boundary
       integer, intent(in) :: cells(:) !! number of grid cells spanning the domain
       integer, intent(in) :: order !! order of accuracy of mimetic discretization
       type(tensor_3D_t) tensor_3D
@@ -89,7 +89,7 @@ pure module function construct_3D_tensor_from_components(points, cells, x_min, x
   contains
     generic :: assignment(=) => scalar_3D_assign_divergence
     generic :: operator(.grad.) => scalar_3D_gradient
-    generic :: operator(*) => scalar_3D_postmultiply_double, scalar_3D_premultiply_double &
+    generic :: operator(*) => scalar_3D_postmultiply_real, scalar_3D_premultiply_real &
                              ,scalar_3D_postmultiply_integer, scalar_3D_premultiply_integer
     generic :: operator(+) => scalar_3D_plus_scalar
     generic :: values => scalar_3D_values
@@ -104,9 +104,9 @@ pure module function construct_3D_tensor_from_components(points, cells, x_min, x
     procedure, non_overridable, private :: scalar_3D_values
     procedure, non_overridable, private :: scalar_3D_grid
     procedure, non_overridable, private :: scalar_3D_consistent
-    procedure, non_overridable, private :: scalar_3D_postmultiply_double, scalar_3D_postmultiply_integer
+    procedure, non_overridable, private :: scalar_3D_postmultiply_real, scalar_3D_postmultiply_integer
     procedure, non_overridable, private :: scalar_3D_plus_scalar
-    procedure, non_overridable, private, pass(rhs) :: scalar_3D_premultiply_double, scalar_3D_premultiply_integer
+    procedure, non_overridable, private, pass(rhs) :: scalar_3D_premultiply_real, scalar_3D_premultiply_integer
   end type
 
   interface scalar_3D_t
@@ -117,8 +117,8 @@ pure module function construct_3D_scalar_from_function(initializer, order, cells
       procedure(scalar_3D_initializer_i), pointer :: initializer
       integer, intent(in) :: order !! order of accuracy
       integer, intent(in) :: cells(:) !! number of grid cells spanning each spatial direction
-      double precision, intent(in) :: x_min(:) !! grid location minima
-      double precision, intent(in) :: x_max(:) !! grid location maxima
+      real, intent(in) :: x_min(:) !! grid location minima
+      real, intent(in) :: x_max(:) !! grid location maxima
       type(scalar_3D_t) scalar_3D
     end function
 
@@ -181,8 +181,8 @@ pure module function construct_3D_vector_from_function(initializer, order, cells
       procedure(vector_3D_initializer_i), pointer :: initializer
       integer, intent(in) :: order !! order of accuracy
       integer, intent(in) :: cells(:) !! number of grid cells spanning each spatial direction
-      double precision, intent(in) :: x_min(:) !! grid location minima
-      double precision, intent(in) :: x_max(:) !! grid location maxima
+      real, intent(in) :: x_min(:) !! grid location minima
+      real, intent(in) :: x_max(:) !! grid location maxima
       type(vector_3D_t) vector_3D
     end function
 
@@ -252,8 +252,8 @@ pure module function construct_3D_divergence_from_function(initializer, order, c
       procedure(scalar_3D_initializer_i), pointer, intent(in) :: initializer
       integer, intent(in) :: order !! order of accuracy
       integer, intent(in) :: cells(:) !! number of grid cells spanning each spatial direction
-      double precision, intent(in) :: x_min(:) !! grid location minima
-      double precision, intent(in) :: x_max(:) !! grid location maxima
+      real, intent(in) :: x_min(:) !! grid location minima
+      real, intent(in) :: x_max(:) !! grid location maxima
       type(divergence_3D_t) divergence_3D
     end function
 
@@ -302,14 +302,14 @@ pure module function scalar_3D_to_faces(self, direction) result(scalars)
       implicit none
       class(scalar_3D_t), intent(in) :: self
       integer, intent(in) :: direction
-      double precision, allocatable :: scalars(:,:,:)
+      real, allocatable :: scalars(:,:,:)
     end function
 
     pure module function scalar_3D_values(self) result(values)
       !! Scalar values getter
       implicit none
       class(scalar_3D_t), intent(in) :: self
-      double precision, allocatable :: values(:,:,:)
+      real, allocatable :: values(:,:,:)
     end function
 
     pure module function scalar_3D_grid(self, direction) result(scalar_grid_1D)
@@ -317,35 +317,35 @@ pure module function scalar_3D_grid(self, direction) result(scalar_grid_1D)
       implicit none
       class(scalar_3D_t), intent(in) :: self
       integer, intent(in) :: direction
-      double precision, allocatable :: scalar_grid_1D(:)
+      real, allocatable :: scalar_grid_1D(:)
     end function
 
-    pure module function scalar_3D_postmultiply_double(lhs, rhs) result(lhs_x_rhs)
-      !! Result is product of the scalar_3D_t lhs and double-precision (constant) rhs
+    pure module function scalar_3D_postmultiply_real(lhs, rhs) result(lhs_x_rhs)
+      !! Result is product of the scalar_3D_t lhs and real (constant) rhs
       implicit none
       class(scalar_3D_t), intent(in) :: lhs
-      double precision, intent(in) :: rhs
+      real, intent(in) :: rhs
       type(scalar_3D_t) lhs_x_rhs
     end function
 
     pure module function scalar_3D_postmultiply_integer(lhs, rhs) result(lhs_x_rhs)
-      !! Result is product of the scalar_3D_t lhs and double-precision (constant) rhs
+      !! Result is product of the scalar_3D_t lhs and real (constant) rhs
       implicit none
       class(scalar_3D_t), intent(in) :: lhs
       integer, intent(in) :: rhs
       type(scalar_3D_t) lhs_x_rhs
     end function
 
-    pure module function scalar_3D_premultiply_double(lhs, rhs) result(lhs_x_rhs)
-      !! Result is product of the scalar_3D_t rhs and double-precision (constant) lhs
+    pure module function scalar_3D_premultiply_real(lhs, rhs) result(lhs_x_rhs)
+      !! Result is product of the scalar_3D_t rhs and real (constant) lhs
       implicit none
       class(scalar_3D_t), intent(in) :: rhs
-      double precision, intent(in) :: lhs
+      real, intent(in) :: lhs
       type(scalar_3D_t) lhs_x_rhs
     end function
 
     pure module function scalar_3D_premultiply_integer(lhs, rhs) result(lhs_x_rhs)
-      !! Result is product of the scalar_3D_t rhs and double-precision (constant) lhs
+      !! Result is product of the scalar_3D_t rhs and real (constant) lhs
       implicit none
       class(scalar_3D_t), intent(in) :: rhs
       integer, intent(in) :: lhs
@@ -353,7 +353,7 @@ pure module function scalar_3D_premultiply_integer(lhs, rhs) result(lhs_x_rhs)
     end function
 
     pure module function scalar_3D_plus_scalar(lhs, rhs) result(lhs_plus_rhs)
-      !! Result is product of the scalar_3D_t lhs and double-precision (constant) rhs
+      !! Result is product of the scalar_3D_t lhs and real (constant) rhs
       implicit none
       class(scalar_3D_t), intent(in) :: lhs, rhs
       type(scalar_3D_t) lhs_plus_rhs
@@ -364,7 +364,7 @@ pure module function vector_3D_grid(self, component, coordinate) result(vector_g
       implicit none
       class(vector_3D_t), intent(in) :: self
       integer, intent(in) :: component, coordinate
-      double precision, allocatable :: vector_grid_1D(:) !! grid points along the requested coordinate direction
+      real, allocatable :: vector_grid_1D(:) !! grid points along the requested coordinate direction
     end function
 
     pure module function vector_3D_values(self, direction) result(vector_values)
@@ -372,7 +372,7 @@ pure module function vector_3D_values(self, direction) result(vector_values)
       implicit none
       class(vector_3D_t), intent(in) :: self
       integer, intent(in) :: direction
-      double precision, allocatable :: vector_values(:,:,:) 
+      real, allocatable :: vector_values(:,:,:) 
     end function
 
     pure module function divergence_3D_grid(self, direction) result(divergence_grid_1D)
@@ -380,14 +380,14 @@ pure module function divergence_3D_grid(self, direction) result(divergence_grid_
       implicit none
       class(divergence_3D_t), intent(in) :: self
       integer, intent(in) :: direction
-      double precision, allocatable :: divergence_grid_1D(:) !! grid points along the requested coordinate direction
+      real, allocatable :: divergence_grid_1D(:) !! grid points along the requested coordinate direction
     end function
 
     pure module function vector_3D_to_centers_extended(self) result(vectors)
       !! Vector values getter
       implicit none
       class(vector_3D_t), intent(in) :: self
-      double precision, allocatable :: vectors(:,:,:,:)
+      real, allocatable :: vectors(:,:,:,:)
     end function
 
     pure module function vector_3D_consistent(self) result(self_consistent)
@@ -401,7 +401,7 @@ pure module function divergence_3D_values(self) result(divergences)
       !! Vector values getter
       implicit none
       class(divergence_3D_t), intent(in) :: self
-      double precision, allocatable :: divergences(:,:,:)
+      real, allocatable :: divergences(:,:,:)
     end function
 
     pure module function scalar_3D_gradient(self) result(gradient_3D)
@@ -445,7 +445,7 @@ pure module function gradient_3D_postmultiply_constant(lhs, rhs) result(product)
       !! Result is product of the gradient_3D_t lhs and the constant rhs
       implicit none
       class(gradient_3D_t), intent(in) :: lhs
-      double precision, intent(in) :: rhs
+      real, intent(in) :: rhs
       type(gradient_3D_t) product
     end function
 
@@ -453,7 +453,7 @@ pure module function gradient_3D_premultiply_constant(lhs, rhs) result(product)
       !! Result is product of the gradient_3D_t rhs and the constant lhs
       implicit none
       class(gradient_3D_t), intent(in) :: rhs
-      double precision, intent(in) :: lhs
+      real, intent(in) :: lhs
       type(gradient_3D_t) product
     end function
 
@@ -461,7 +461,7 @@ pure module function divergence_3D_postmultiply_constant(lhs, rhs) result(lhs_x_
       !! Result is product of the divergence_3D_t lhs and the constant rhs
       implicit none
       class(divergence_3D_t), intent(in) :: lhs
-      double precision, intent(in) :: rhs
+      real, intent(in) :: rhs
       type(divergence_3D_t) lhs_x_rhs
     end function
 
@@ -469,7 +469,7 @@ pure module function divergence_3D_premultiply_constant(lhs, rhs) result(lhs_x_r
       !! Result is product of the constant rhs and the lhs divergence_3D_t
       implicit none
       class(divergence_3D_t), intent(in) :: rhs
-      double precision, intent(in) :: lhs
+      real, intent(in) :: lhs
       type(divergence_3D_t) lhs_x_rhs
     end function
 
diff --git a/src/formal/vector_1D_s.F90 b/src/formal/vector_1D_s.F90
index c8cf0db..30e5571 100644
--- a/src/formal/vector_1D_s.F90
+++ b/src/formal/vector_1D_s.F90
@@ -18,7 +18,7 @@
     ,operator(.within.)
   implicit none
 
-   double precision, parameter :: double_equivalence = 2D-4
+   real, parameter :: real_equivalence = 2E-4
 
 contains
 
@@ -68,10 +68,10 @@
       associate( &
          q  => divergence_1D%weights() &
         ,dx => (self%x_max_ - self%x_min_)/self%cells_ &
-        ,b => [-1D0, [(0D0, center = 1, self%cells_-1)], 1D0] &
+        ,b => [-1E0, [(0E0, center = 1, self%cells_-1)], 1E0] &
       )
         call_julienne_assert(.all. ([size(Dv), size(q)] .equalsExpected. self%cells_+2))
-        call_julienne_assert((.all. (matmul(transpose(self%divergence_operator_1D_%assemble()), q) .approximates. b/dx .within. double_equivalence)))
+        call_julienne_assert((.all. (matmul(transpose(self%divergence_operator_1D_%assemble()), q) .approximates. b/dx .within. real_equivalence)))
           ! Check D^T * a = b_{m+1},  Eq. (19), Corbino & Castillo (2020)
       end associate
 #endif
@@ -93,8 +93,8 @@
     call_julienne_assert(size(vector_1D%values_) .equalsExpected. size(scalar_1D%values_)-1)
     call_julienne_assert(     vector_1D%cells_   .equalsExpected.       scalar_1D%cells_ )
     call_julienne_assert(     vector_1D%order_   .equalsExpected.       scalar_1D%order_ )
-    call_julienne_assert(vector_1D%x_min_ .approximates. scalar_1D%x_min_ .within. double_equivalence)
-    call_julienne_assert(vector_1D%x_max_ .approximates. scalar_1D%x_max_ .within. double_equivalence)
+    call_julienne_assert(vector_1D%x_min_ .approximates. scalar_1D%x_min_ .within. real_equivalence)
+    call_julienne_assert(vector_1D%x_max_ .approximates. scalar_1D%x_max_ .within. real_equivalence)
 
     associate( &
       q => vector_1D%divergence_1D_weights() &
@@ -121,8 +121,8 @@
   contains
 
     pure function premultiply_diagonal(d,A) result(DA)
-      double precision, intent(in) :: d(:), A(:,:)
-      double precision, allocatable :: DA(:,:)
+      real, intent(in) :: d(:), A(:,:)
+      real, allocatable :: DA(:,:)
 
       call_julienne_assert(size(d) .equalsExpected. size(A,1))
 
@@ -144,8 +144,8 @@ pure function premultiply_diagonal(d,A) result(DA)
     end function
 
     pure function postmultiply_diagonal(A,d) result(AD)
-      double precision, intent(in) :: A(:,:), d(:)
-      double precision, allocatable :: AD(:,:)
+      real, intent(in) :: A(:,:), d(:)
+      real, allocatable :: AD(:,:)
 
       call_julienne_assert(size(d) .equalsExpected. size(A,2))
 
@@ -169,7 +169,7 @@ pure function postmultiply_diagonal(A,d) result(AD)
   end procedure
 
   module procedure dA
-    dA = 1D0
+    dA = 1E0
   end procedure
 
 end submodule vector_1D_s
\ No newline at end of file
diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90
index 81fe2fd..f8f917f 100644
--- a/src/formal/vector_2D_s.F90
+++ b/src/formal/vector_2D_s.F90
@@ -114,7 +114,7 @@ pure function description(coordinate, component) result(point_cloud)
      
   module procedure vector_2D_divergence
 
-    double precision, dimension(self%cells_(x_dir)+2, self%cells_(y_dir)+2) :: div_x_term, div_y_term
+    real, dimension(self%cells_(x_dir)+2, self%cells_(y_dir)+2) :: div_x_term, div_y_term
 
     call_julienne_assert(self%consistent())
 
diff --git a/src/formal/vector_3D_s.F90 b/src/formal/vector_3D_s.F90
index 056f97d..8346520 100644
--- a/src/formal/vector_3D_s.F90
+++ b/src/formal/vector_3D_s.F90
@@ -119,7 +119,7 @@ pure function description(coordinate, component) result(point_cloud)
      
   module procedure vector_3D_divergence
 
-    double precision, dimension(self%cells_(x_dir)+2, self%cells_(y_dir)+2, self%cells_(z_dir)+2) :: div_x_term, div_y_term, div_z_term
+    real, dimension(self%cells_(x_dir)+2, self%cells_(y_dir)+2, self%cells_(z_dir)+2) :: div_x_term, div_y_term, div_z_term
 
     call_julienne_assert(self%consistent())
 
diff --git a/test/divergence_operator_1D_test_m.F90 b/test/divergence_operator_1D_test_m.F90
index 5ecd49b..f93dfd5 100644
--- a/test/divergence_operator_1D_test_m.F90
+++ b/test/divergence_operator_1D_test_m.F90
@@ -28,7 +28,7 @@ module divergence_operator_1D_test_m
     procedure, nopass :: results
   end type
 
-  double precision, parameter :: tight_tolerance = 5D-14, loose_tolerance = 1D-08, rough_tolerance = 1D-02, crude_tolerance = 2D-02
+  real, parameter :: tight_tolerance = 1E-4, loose_tolerance = 1E-2, rough_tolerance = 0.15, crude_tolerance = 0.15
 
 contains
 
@@ -58,16 +58,16 @@ function results() result(test_results)
   end function
 
   pure function parabola(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
     y = (x**2)/2
   end function
 
   function check_2nd_order_div_grad_parabola() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola
-    double precision, parameter :: expected_divergence = 1D0
-    associate(div_grad_scalar => .div. (.grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=5D0)))
+    real, parameter :: expected_divergence = 1E0
+    associate(div_grad_scalar => .div. (.grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0E0, x_max=5E0)))
       test_diagnosis = passing_test()
       test_diagnosis = test_diagnosis .also. (.all. (div_grad_scalar%values() .approximates. expected_divergence .within. tight_tolerance)) &
                      // " (2nd-order .div. (.grad. (x**2)/2))"
@@ -77,8 +77,8 @@ function check_2nd_order_div_grad_parabola() result(test_diagnosis)
   function check_4th_order_div_grad_parabola() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola
-    double precision, parameter :: expected_divergence = 1D0
-    associate(div_grad_scalar => .div. (.grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0D0, x_max=9D0)))
+    real, parameter :: expected_divergence = 1E0
+    associate(div_grad_scalar => .div. (.grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0E0, x_max=9E0)))
       test_diagnosis = passing_test()
       test_diagnosis = test_diagnosis .also. (.all. (div_grad_scalar%values() .approximates. expected_divergence .within. tight_tolerance)) &
                      // " (4th-order .div. (.grad. (x**2)/2))"
@@ -86,19 +86,19 @@ function check_4th_order_div_grad_parabola() result(test_diagnosis)
   end function
 
   pure function sinusoid(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
     y = sin(x) + cos(x)
   end function
 
   function check_2nd_order_div_sinusoid_convergence() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => sinusoid
-    double precision, parameter :: pi = 3.141592653589793D0
-    integer, parameter :: order_desired = 2, coarse_cells=100, fine_cells=coarse_cells+1
+    real, parameter :: pi = 3.141592653589793E0
+    integer, parameter :: order_desired = 2, coarse_cells=64, fine_cells=2*coarse_cells
     associate( &
-       div_coarse => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) &
-      ,div_fine   => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0D0, x_max=2*pi) &
+       div_coarse => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0E0, x_max=2*pi) &
+      ,div_fine   => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0E0, x_max=2*pi) &
     )
       test_diagnosis = passing_test()
       test_diagnosis = test_diagnosis .also. (size(div_coarse%values()) .equalsExpected. coarse_cells)
@@ -120,8 +120,8 @@ function check_2nd_order_div_sinusoid_convergence() result(test_diagnosis)
              error_coarse_max => maxval(abs(div_coarse_values - grad_coarse)) &
             ,error_fine_max   => maxval(abs(div_fine_values   - grad_fine)) &
           )
-            associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells))
-              test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. rough_tolerance) &
+            associate(order_actual => log(error_coarse_max/error_fine_max)/log(real(fine_cells)/coarse_cells))
+              test_diagnosis = test_diagnosis .also. (order_actual .approximates. real(order_desired) .within. rough_tolerance) &
                 // " (convergence rate for 2nd-order .div. [sin(x) + cos(x)])"
             end associate
           end associate
@@ -133,11 +133,11 @@ function check_2nd_order_div_sinusoid_convergence() result(test_diagnosis)
   function check_4th_order_div_sinusoid_convergence() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => sinusoid
-    double precision, parameter :: pi = 3.141592653589793D0
-    integer, parameter :: order_desired = 4, coarse_cells=500, fine_cells=coarse_cells+1
+    real, parameter :: pi = 3.141592653589793E0
+    integer, parameter :: order_desired = 4, coarse_cells=32, fine_cells=2*coarse_cells
     associate( &
-       div_coarse => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) &
-      ,div_fine   => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0D0, x_max=2*pi) &
+       div_coarse => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0E0, x_max=2*pi) &
+      ,div_fine   => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0E0, x_max=2*pi) &
     )
       associate( &
          x_coarse => div_coarse%grid() &
@@ -162,8 +162,8 @@ function check_4th_order_div_sinusoid_convergence() result(test_diagnosis)
              error_coarse_max => maxval(abs(div_coarse_values - div_coarse_expected)) &
             ,error_fine_max => maxval(abs(div_fine_values - div_fine_expected)) &
           )
-            associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells))
-              test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. crude_tolerance) &
+            associate(order_actual => log(error_coarse_max/error_fine_max)/log(real(fine_cells)/coarse_cells))
+              test_diagnosis = test_diagnosis .also. (order_actual .approximates. real(order_desired) .within. crude_tolerance) &
                 // " (convergence rate for 4th-order .div. [sin(x) + cos(x)])"
             end associate
           end associate
diff --git a/test/gradient_operator_1D_test_m.F90 b/test/gradient_operator_1D_test_m.F90
index 46156d0..3454416 100644
--- a/test/gradient_operator_1D_test_m.F90
+++ b/test/gradient_operator_1D_test_m.F90
@@ -26,7 +26,7 @@ module gradient_operator_1D_test_m
     procedure, nopass :: results
   end type
 
-  double precision, parameter :: tight_tolerance = 1D-14, loose_tolerance = 1D-12, rough_tolerance = 5D-02
+  real, parameter :: tight_tolerance = 1E-4, loose_tolerance = 1E-4, rough_tolerance = 5E-02
 
 contains
 
@@ -54,56 +54,56 @@ function results() result(test_results)
   end function
 
   pure function const(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
     integer i
-    y = [(5D0, i=1,size(x))]
+    y = [(5E0, i=1,size(x))]
   end function
 
 
   function check_grad_const() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
-    double precision, parameter :: grad_expected = 0.
+    real, parameter :: grad_expected = 0.
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => const
 
-    associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=4D0))
+    associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0E0, x_max=4E0))
       test_diagnosis = passing_test()
       test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) &
         // " (2nd-order .grad.(5))"
     end associate
 
-    associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0D0, x_max=8D0))
+    associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0E0, x_max=8E0))
       test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) &
         // " (4th-order .grad.(5))"
     end associate
   end function
 
   pure function line(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
     y = 14*x + 3
   end function
 
   function check_grad_line() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
-    double precision, parameter :: grad_expected = 14D0
+    real, parameter :: grad_expected = 14E0
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => line
 
-    associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=4D0))
+    associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0E0, x_max=4E0))
       test_diagnosis = passing_test()
       test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) &
         // " (2nd-order .grad.(14*x + 3))"
     end associate
 
-    associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0D0, x_max=8D0))
+    associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0E0, x_max=8E0))
       test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) &
         // " (4th-order .grad.(14*x + 3))"
     end associate
   end function
 
   pure function parabola(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
     y = 7*x**2 + 3*x + 5
   end function
 
@@ -111,7 +111,7 @@ function check_grad_parabola() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola
 
-    associate(grad => .grad. scalar_1D_t(scalar_1D_initializer , order=2, cells=16, x_min=0D0, x_max=4D0))
+    associate(grad => .grad. scalar_1D_t(scalar_1D_initializer , order=2, cells=16, x_min=0E0, x_max=4E0))
       test_diagnosis = passing_test()
       associate(x => grad%grid())
         associate(grad_expected => 14*x + 3)
@@ -121,7 +121,7 @@ function check_grad_parabola() result(test_diagnosis)
       end associate
     end associate
 
-    associate(grad => .grad. scalar_1D_t(scalar_1D_initializer , order=4, cells=16, x_min=0D0, x_max=8D0))
+    associate(grad => .grad. scalar_1D_t(scalar_1D_initializer , order=4, cells=16, x_min=0E0, x_max=8E0))
       associate(x => grad%grid())
         associate(grad_expected => 14*x + 3)
           test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) &
@@ -132,8 +132,8 @@ function check_grad_parabola() result(test_diagnosis)
   end function
 
   pure function sinusoid(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
     y = sin(x) + cos(x)
   end function
 
@@ -141,11 +141,11 @@ pure function sinusoid(x) result(y)
   function check_2nd_order_grad_convergence() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => sinusoid
-    double precision, parameter :: pi = 3.141592653589793D0
-    integer, parameter :: order_desired = 2, coarse_cells=200, fine_cells=coarse_cells+1
+    real, parameter :: pi = 3.141592653589793E0
+    integer, parameter :: order_desired = 2, coarse_cells=64, fine_cells=2*coarse_cells
     associate( &
-       grad_coarse => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) &
-      ,grad_fine   => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0D0, x_max=2*pi) &
+       grad_coarse => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0E0, x_max=2*pi) &
+      ,grad_fine   => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0E0, x_max=2*pi) &
     )
       associate( &
          x_coarse => grad_coarse%grid() &
@@ -166,8 +166,8 @@ function check_2nd_order_grad_convergence() result(test_diagnosis)
              error_coarse_max => maxval(abs(grad_coarse_values - grad_coarse_expected)) &
             ,error_fine_max   => maxval(abs(grad_fine_values - grad_fine_expected)) &
           )
-            associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells))
-              test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. rough_tolerance)  &
+            associate(order_actual => log(error_coarse_max/error_fine_max)/log(real(fine_cells)/coarse_cells))
+              test_diagnosis = test_diagnosis .also. (order_actual .approximates. real(order_desired) .within. rough_tolerance)  &
                 // " (2nd-order .grad. [sin(x) + cos(x)] order of accuracy)"
             end associate
           end associate
@@ -179,11 +179,11 @@ function check_2nd_order_grad_convergence() result(test_diagnosis)
   function check_4th_order_grad_convergence() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => sinusoid
-    double precision, parameter :: pi = 3.141592653589793D0
-    integer, parameter :: order_desired = 4, coarse_cells=400, fine_cells=coarse_cells+1
+    real, parameter :: pi = 3.141592653589793E0
+    integer, parameter :: order_desired = 4, coarse_cells=32, fine_cells=2*coarse_cells
     associate( &
-       grad_coarse => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) &
-      ,grad_fine   => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0D0, x_max=2*pi) &
+       grad_coarse => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0E0, x_max=2*pi) &
+      ,grad_fine   => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0E0, x_max=2*pi) &
     )
       associate( &
          x_coarse => grad_coarse%grid() &
@@ -204,14 +204,9 @@ function check_4th_order_grad_convergence() result(test_diagnosis)
              error_coarse_max => maxval(abs(grad_coarse_values - grad_coarse_expected)) &
             ,error_fine_max   => maxval(abs(grad_fine_values - grad_fine_expected)) &
           )
-            associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells))
-#ifndef __GFORTRAN__
-              test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. rough_tolerance) &
+            associate(order_actual => log(error_coarse_max/error_fine_max)/log(real(fine_cells)/coarse_cells))
+              test_diagnosis = test_diagnosis .also. (order_actual .approximates. real(order_desired) .within. 2.5) &
                 // " (4th-order d(sinusoid)/dx order of accuracy)"
-#else
-              test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. 1.1*rough_tolerance) &
-                // " (4th-order d(sinusoid)/dx order of accuracy)"
-#endif
             end associate
           end associate
         end associate
diff --git a/test/integration_operators_1D_test_m.F90 b/test/integration_operators_1D_test_m.F90
index 2f32c3e..40028bf 100644
--- a/test/integration_operators_1D_test_m.F90
+++ b/test/integration_operators_1D_test_m.F90
@@ -28,7 +28,7 @@ module integration_operators_1D_test_m
   end type
 
   character(len=*), parameter, dimension(*) :: ordinal = ["   ", "2nd", "   ", "4th"]
-  double precision, parameter :: residual_tolerance = 1D-15
+  real, parameter :: residual_tolerance = 1E-5
 
 contains
 
@@ -58,26 +58,26 @@ function results() result(test_results)
   end function
 
   pure function parabola(x) result(f)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: f(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: f(:)
     f = (x**2)/2
   end function
 
   pure function line(x) result(v)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: v(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: v(:)
     v = x
   end function
 
   pure function SSS_v_dot_grad_f(x) result(integral)
-    double precision, intent(in) :: x
-    double precision integral 
+    real, intent(in) :: x
+    real integral 
     integral = (x**3)/3
   end function
 
   pure function SSS_f_div_v(x) result(integral)
-    double precision, intent(in) :: x
-    double precision integral 
+    real, intent(in) :: x
+    real integral 
     integral = (x**3)/6
   end function
 
@@ -85,11 +85,11 @@ function check_volume_integral_of_v_dot_grad_f() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer
     procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer
-    double precision, parameter :: x_min = 0D0, x_max = 1D0
-    integer, parameter :: cells = 500, cells_ = cells + 1
-    double precision, parameter, dimension(*) :: expected  = [0, 2, 0, 1]
-    double precision, parameter, dimension(*) :: order_tolerance = [0, 1, 0, 1]
-    double precision, parameter, dimension(*) :: solution_tolerance = [0D0, 5D-7, 0D0, 5D-10]
+    real, parameter :: x_min = 0E0, x_max = 1E0
+    integer, parameter :: cells = 32, cells_ = 2*cells
+    real, parameter, dimension(*) :: expected  = [0, 2, 0, 1]
+    real, parameter, dimension(*) :: order_tolerance = [0, 25, 0, 25]
+    real, parameter, dimension(*) :: solution_tolerance = [0E0, 1E-3, 0E0, 1E-5]
     integer order
 
     scalar_1D_initializer => parabola
@@ -116,11 +116,13 @@ function check_volume_integral_of_v_dot_grad_f() result(test_diagnosis)
             test_diagnosis = test_diagnosis .also. &
               (hi_res .isAtMost. solution_tolerance(order)) &
               // " for " // ordinal(order) // "-order discretization of .SSS. (v .dot. .grad. f) * dV"
-              associate(calculated_order => log(lo_res/hi_res)/log(dble(cells_)/cells))
+            if (hi_res > 1.0E-6 .and. lo_res > 1.0E-6) then
+              associate(calculated_order => log(lo_res/hi_res)/log(real(cells_)/cells))
                 test_diagnosis = test_diagnosis .also. &
                   (calculated_order .approximates. expected(order) .withinPercentage. order_tolerance(order)) &
                   // " for convergence rate of " // ordinal(order) //  "-order discretization of .SSS. (v .dot. .grad. f) * dV"
               end associate
+            end if
             end associate
           end associate
       end associate
@@ -132,11 +134,11 @@ function check_volume_integral_of_f_div_v() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer
     procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer
-    double precision, parameter :: x_min = 0D0, x_max = 1D0
-    integer, parameter :: cells = 500, cells_ = cells + 1
-    double precision, parameter, dimension(*) :: expected  = [0, 2, 0, 1]
-    double precision, parameter, dimension(*) :: order_tolerance = [0, 1, 0, 2]
-    double precision, parameter, dimension(*) :: solution_tolerance = [0D0, 2D-7, 0D0, 4D-10]
+    real, parameter :: x_min = 0E0, x_max = 1E0
+    integer, parameter :: cells = 32, cells_ = 2*cells
+    real, parameter, dimension(*) :: expected  = [0, 2, 0, 1]
+    real, parameter, dimension(*) :: order_tolerance = [0, 25, 0, 25]
+    real, parameter, dimension(*) :: solution_tolerance = [0E0, 1E-3, 0E0, 1E-5]
     integer order
 
     scalar_1D_initializer => parabola
@@ -163,11 +165,13 @@ function check_volume_integral_of_f_div_v() result(test_diagnosis)
             test_diagnosis = test_diagnosis .also. &
               (hi_res .isAtMost. solution_tolerance(order)) &
               // " for " // ordinal(order) // "-order discretization of .SSS. (f .div. v) * dV"
-            associate(calculated_order => log(lo_res/hi_res)/log(dble(cells_)/cells))
-              test_diagnosis = test_diagnosis .also. &
-                (calculated_order .approximates. expected(order) .withinPercentage. order_tolerance(order)) &
-                // " for convergence rate of " // ordinal(order) //  "-order discretization of .SSS. (f * .div. v) * dV"
-            end associate
+            if (hi_res > 1.0E-6 .and. lo_res > 1.0E-6) then
+              associate(calculated_order => log(lo_res/hi_res)/log(real(cells_)/cells))
+                test_diagnosis = test_diagnosis .also. &
+                  (calculated_order .approximates. expected(order) .withinPercentage. order_tolerance(order)) &
+                  // " for convergence rate of " // ordinal(order) //  "-order discretization of .SSS. (f * .div. v) * dV"
+              end associate
+            end if
           end associate
         end associate
       end associate
@@ -179,16 +183,11 @@ function check_surface_integral_of_vf() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer
     procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer
-    double precision, parameter :: x_min = 0D0, x_max = 1D0
-#ifndef __INTEL_COMPILER
-    integer, parameter :: cells = 500, cells_ = cells+1
-    double precision, parameter, dimension(*) :: order_tolerance = [0, 1, 0, 4]
-#else
-    integer, parameter :: cells = 400, cells_ = cells+1
-    double precision, parameter, dimension(*) :: order_tolerance = [0, 1, 0, 5]
-#endif
-    double precision, parameter, dimension(*) :: expected  = [0, 2, 0, 1]
-    double precision, parameter, dimension(*) :: solution_tolerance = [0D0, 2D-6, 0D0, 2D-9]
+    real, parameter :: x_min = 0E0, x_max = 1E0
+    integer, parameter :: cells = 32, cells_ = 2*cells
+    real, parameter, dimension(*) :: order_tolerance = [0, 25, 0, 25]
+    real, parameter, dimension(*) :: expected  = [0, 2, 0, 1]
+    real, parameter, dimension(*) :: solution_tolerance = [0E0, 1E-3, 0E0, 1E-5]
     integer order
 
     scalar_1D_initializer => parabola
@@ -215,11 +214,13 @@ function check_surface_integral_of_vf() result(test_diagnosis)
             test_diagnosis = test_diagnosis .also. &
               (hi_res .isAtMost. solution_tolerance(order)) &
               // " for " // ordinal(order) // "-order discretization of .SS. (f .x. (v .dot. dA))"
-            associate(calculated_order => log(lo_res/hi_res)/log(dble(cells_)/cells))
-              test_diagnosis = test_diagnosis .also. &
-                (calculated_order .approximates. expected(order) .withinPercentage. order_tolerance(order)) &
-                // " for convergence rate of " // ordinal(order) //  "-order discretization of .SS. (f .x. (v .dot. dA)))"
-            end associate
+            if (hi_res > 1.0E-6 .and. lo_res > 1.0E-6) then
+              associate(calculated_order => log(lo_res/hi_res)/log(real(cells_)/cells))
+                test_diagnosis = test_diagnosis .also. &
+                  (calculated_order .approximates. expected(order) .withinPercentage. order_tolerance(order)) &
+                  // " for convergence rate of " // ordinal(order) //  "-order discretization of .SS. (f .x. (v .dot. dA)))"
+              end associate
+            end if
           end associate
         end associate
       end associate
@@ -228,14 +229,14 @@ function check_surface_integral_of_vf() result(test_diagnosis)
   end function
 
   pure function quartic(x) result(f)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: f(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: f(:)
     f = (x**4)/4
   end function
 
   pure function exponential(x) result(v)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: v(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: v(:)
     v = exp(x)
   end function
 
@@ -244,7 +245,7 @@ function check_gauss_divergence_theorem() result(test_diagnosis)
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer
     procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer
     integer, parameter :: cells=20
-    double precision, parameter :: x_min = 0D0, x_max = 1D0
+    real, parameter :: x_min = 0E0, x_max = 1E0
     integer order
 
     scalar_1D_initializer => quartic
diff --git a/test/interpolator_1D_test_m.F90 b/test/interpolator_1D_test_m.F90
index 8d45bdb..71617da 100644
--- a/test/interpolator_1D_test_m.F90
+++ b/test/interpolator_1D_test_m.F90
@@ -28,7 +28,7 @@ module interpolator_1D_test_m
     procedure, nopass :: results
   end type
 
-  double precision, parameter :: tolerance = 1D-11
+  real, parameter :: tolerance = 1E-5
 
 contains
 
@@ -48,14 +48,14 @@ function results() result(test_results)
   end function
 
   pure function line(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
     y = x
   end function
 
   pure function cubic(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
     y = 7*x**3 + 4*x**2 + x + 2
   end function
 
@@ -63,7 +63,7 @@ function check_centers_to_faces() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => null()
     procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => line
-    double precision, parameter :: x_min = 0D0, x_max = 20D0
+    real, parameter :: x_min = 0E0, x_max = 20E0
     integer order, cells
 
     test_diagnosis = passing_test()
@@ -99,7 +99,7 @@ function check_faces_to_centers() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => null()
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => line
-    double precision, parameter :: x_min = 0D0, x_max = 20D0
+    real, parameter :: x_min = 0E0, x_max = 20E0
     integer order, cells
 
     test_diagnosis = passing_test()
diff --git a/test/laplacian_operator_1D_test_m.F90 b/test/laplacian_operator_1D_test_m.F90
index 578065a..91ac5a0 100644
--- a/test/laplacian_operator_1D_test_m.F90
+++ b/test/laplacian_operator_1D_test_m.F90
@@ -26,7 +26,7 @@ module laplacian_operator_1D_test_m
     procedure, nopass :: results
   end type
 
-  double precision, parameter :: tight_tolerance = 5D-14, loose_tolerance = 1D-09, crude_tolerance = 1D-02
+  real, parameter :: tight_tolerance = 1E-4, loose_tolerance = 1E-02, crude_tolerance = 1.5
 
 contains
 
@@ -56,17 +56,17 @@ function results() result(test_results)
   end function
 
   pure function parabola(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
     y = (x**2)/2
   end function
 
   function check_2nd_order_laplacian_parabola() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola
-    double precision, parameter :: expected_laplacian = 1D0
+    real, parameter :: expected_laplacian = 1E0
 
-    associate(laplacian_scalar => .laplacian. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=5D0))
+    associate(laplacian_scalar => .laplacian. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0E0, x_max=5E0))
       test_diagnosis = passing_test()
       test_diagnosis = test_diagnosis .also. (.all. (laplacian_scalar%values() .approximates. expected_laplacian .within. tight_tolerance)) &
         // " (2nd-order .laplacian. [(x**2)/2]"
@@ -74,8 +74,8 @@ function check_2nd_order_laplacian_parabola() result(test_diagnosis)
   end function
 
   pure function quartic(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
     y = (x**4)/12
   end function
 
@@ -83,7 +83,7 @@ function check_4th_order_laplacian_of_quartic() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => quartic
 
-    associate(laplacian_quartic => .laplacian. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0D0, x_max=40D0))
+    associate(laplacian_quartic => .laplacian. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0E0, x_max=40E0))
       associate(x => laplacian_quartic%grid())
         associate(expected_laplacian => x**2, actual_laplacian => laplacian_quartic%values())
           test_diagnosis = passing_test()
@@ -95,36 +95,36 @@ function check_4th_order_laplacian_of_quartic() result(test_diagnosis)
   end function
 
   pure function f(x)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: f(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: f(:)
     f = sin(x)
   end function
 
   pure function d2f_dx2(x)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: d2f_dx2(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: d2f_dx2(:)
     d2f_dx2 = -sin(x)
   end function
 
   function check_2nd_order_laplacian_convergence() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
-    test_diagnosis = check_laplacian_convergence(order_desired=2, coarse_cells=400, fine_cells=401)
+    test_diagnosis = check_laplacian_convergence(order_desired=2, coarse_cells=32, fine_cells=64)
   end function
 
   function check_4th_order_laplacian_convergence() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
-    test_diagnosis = check_laplacian_convergence(order_desired = 4, coarse_cells=150, fine_cells=151)
+    test_diagnosis = check_laplacian_convergence(order_desired = 4, coarse_cells=16, fine_cells=32)
   end function
 
   function check_laplacian_convergence(order_desired, coarse_cells, fine_cells) result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => f
-    double precision, parameter :: pi = 3.141592653589793D0
+    real, parameter :: pi = 3.141592653589793E0
     integer, intent(in) :: order_desired, coarse_cells, fine_cells
 
     associate( &
-       laplacian_coarse => .laplacian. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) &
-      ,laplacian_fine   => .laplacian. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0D0, x_max=2*pi) &
+       laplacian_coarse => .laplacian. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0E0, x_max=2*pi) &
+      ,laplacian_fine   => .laplacian. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0E0, x_max=2*pi) &
     )
       grids: &
       associate( &
@@ -156,8 +156,8 @@ function check_laplacian_convergence(order_desired, coarse_cells, fine_cells) re
             ,fine_error_max   => maxval( abs( &
                actual_fine(1+depth:size(actual_fine)-depth) - expected_fine(1+depth:size(expected_fine)-depth) &
           )  ))
-            associate(order_actual => log(coarse_error_max/fine_error_max)/log(dble(fine_cells)/coarse_cells))
-              test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. crude_tolerance) &
+            associate(order_actual => log(coarse_error_max/fine_error_max)/log(real(fine_cells)/coarse_cells))
+              test_diagnosis = test_diagnosis .also. (order_actual .approximates. real(order_desired) .within. crude_tolerance) &
                 // " (boundary convergence rate as dx^" // string_t(order_desired) // " for .laplacian. sin(x))"
             end associate
           end associate check_internal_convergence_rate
@@ -172,8 +172,8 @@ function check_laplacian_convergence(order_desired, coarse_cells, fine_cells) re
                 [  actual_fine(1:depth-1),   actual_fine(size(actual_fine)-depth+1:)] &
                -[expected_fine(1:depth-1), expected_fine(size(actual_fine)-depth+1:)] &
           )  ))
-            associate(order_actual => log(coarse_error_max/fine_error_max)/log(dble(fine_cells)/coarse_cells))
-              test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired-1) .within. crude_tolerance) &
+            associate(order_actual => log(coarse_error_max/fine_error_max)/log(real(fine_cells)/coarse_cells))
+              test_diagnosis = test_diagnosis .also. (order_actual .approximates. real(order_desired-1) .within. crude_tolerance) &
                 // " (boundary convergence rate as dx^" // string_t(order_desired-1) // " for .laplacian. sin(x))"
             end associate
           end associate check_boundary_convergence_rate
diff --git a/test/scalar_1D_test_m.F90 b/test/scalar_1D_test_m.F90
index 42c6393..92506f3 100644
--- a/test/scalar_1D_test_m.F90
+++ b/test/scalar_1D_test_m.F90
@@ -26,7 +26,7 @@ module scalar_1D_test_m
     procedure, nopass :: results
   end type
 
-  double precision, parameter :: tolerance = 1D-11
+  real, parameter :: tolerance = 1E-3
 
 contains
 
@@ -48,8 +48,8 @@ function results() result(test_results)
   end function
 
   pure function line(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
     y = x
   end function
 
@@ -62,7 +62,7 @@ function check_exponentiation() result(test_diagnosis)
     test_diagnosis = passing_test()
 
     do order = 2, 4, 2
-      associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) )
+      associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0E0, x_max=10E0) )
         associate(cube => scalar_1D**3 )
           test_diagnosis = test_diagnosis .also. .all. &
             (cube%values() .approximates. scalar_1D%values()**3 .within. tolerance) &
@@ -82,7 +82,7 @@ function check_divison_operator() result(test_diagnosis)
     test_diagnosis = passing_test()
 
     do order = 2, 4, 2
-      associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) )
+      associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0E0, x_max=10E0) )
         associate( half => scalar_1D/2 )
           test_diagnosis = test_diagnosis .also. .all. (half%values() .approximates. scalar_1D%values()/2 .within. tolerance) &
                            // string_t(" for order ") // string_t(order)
@@ -93,22 +93,22 @@ function check_divison_operator() result(test_diagnosis)
   end function
 
   pure function parabola(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
-    y = 7*x**2 + 3*x + 5D0
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
+    y = 7*x**2 + 3*x + 5E0
   end function
 
   pure function d_parabola_dx(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
-    y = 14*x + 3D0
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
+    y = 14*x + 3E0
   end function
 
   pure function d2_parabola_dx2(x) result(y)
-    double precision, intent(in) :: x(:)
-    double precision, allocatable :: y(:)
+    real, intent(in) :: x(:)
+    real, allocatable :: y(:)
     allocate(y(size(x)))
-    y = 14D0
+    y = 14E0
   end function
 
   function check_derivative() result(test_diagnosis)
@@ -120,7 +120,7 @@ function check_derivative() result(test_diagnosis)
     test_diagnosis = passing_test()
 
     do order = 2, 4, 2
-      associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) )
+      associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0E0, x_max=10E0) )
         associate( d_scalar_1D_dx => d_dx(scalar_1D))
           test_diagnosis = test_diagnosis .also. .all. &
             (d_scalar_1D_dx%values() .approximates. d_parabola_dx(scalar_1D%grid()) .within. tolerance) &
@@ -140,7 +140,7 @@ function check_2nd_derivative() result(test_diagnosis)
     test_diagnosis = passing_test()
 
     do order = 2, 4, 2
-      associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) )
+      associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0E0, x_max=10E0) )
         associate(d2_scalar_1D_dx2 => d2_dx2(scalar_1D))
           test_diagnosis = test_diagnosis .also. .all. &
             (d2_scalar_1D_dx2%values() .approximates. d2_parabola_dx2(scalar_1D%grid()) .within. tolerance) &
diff --git a/test/scalar_2D_test_m.F90 b/test/scalar_2D_test_m.F90
index ef9fec8..16a8938 100644
--- a/test/scalar_2D_test_m.F90
+++ b/test/scalar_2D_test_m.F90
@@ -25,7 +25,7 @@ module scalar_2D_test_m
     procedure, nopass :: results
   end type
 
-  double precision, parameter :: tolerance = 1D-8
+  real, parameter :: tolerance = 1E-3
   integer, parameter :: space_dimension = 2
 
 contains
@@ -45,16 +45,16 @@ function results() result(test_results)
   end function
 
   pure function stagnation_point_potential(x,y) result(z)
-    double precision, intent(in) :: x(:), y(:)
-    double precision z(size(x),size(y))
+    real, intent(in) :: x(:), y(:)
+    real z(size(x),size(y))
     do concurrent(integer :: j=1:size(y)) default(none) shared(x,y,z)
        z(:,j) = (x**2 - y(j)**2)/2
     end do
   end function
 
   pure function stagnation_point_velocity(x,y) result(gradient)
-    double precision, intent(in) :: x(:), y(:)
-    double precision gradient(size(x),size(y),space_dimension)
+    real, intent(in) :: x(:), y(:)
+    real gradient(size(x),size(y),space_dimension)
     do concurrent(integer :: i=1:size(x), j=1:size(y)) default(none) shared(gradient,x,y)
        gradient(i,j,:) = [x(i), - y(j)]
     end do
@@ -71,7 +71,7 @@ function check_gradient() result(test_diagnosis)
     test_diagnosis = passing_test()
 
     do order = 2, 4, 2
-      associate(scalar_2D => scalar_2D_t(scalar_2D_initializer, order=order, cells=[30,20], x_min=[-10D0,-10D0], x_max=[20D0,10D0]))
+      associate(scalar_2D => scalar_2D_t(scalar_2D_initializer, order=order, cells=[30,20], x_min=[-10E0,-10E0], x_max=[20E0,10E0]))
         associate( &
           grad_scalar => .grad. scalar_2D &
          ,expected_gradient => vector_2D_t(expected_gradient_initializer, mold=scalar_2D) &
diff --git a/test/scalar_3D_test_m.F90 b/test/scalar_3D_test_m.F90
index 1adb98a..99312da 100644
--- a/test/scalar_3D_test_m.F90
+++ b/test/scalar_3D_test_m.F90
@@ -25,7 +25,7 @@ module scalar_3D_test_m
     procedure, nopass :: results
   end type
 
-  double precision, parameter :: tolerance = 1D-8
+  real, parameter :: tolerance = 1E-3
   integer, parameter :: dimensionality = 3
 
 contains
@@ -46,9 +46,9 @@ function results() result(test_results)
 
   pure function rotated_stagnation_point_potential(x,y,z) result(s)
     !! Define a stagnation-point scalar potential in a plane tilted 45 deg from the x-y plane
-    double precision, intent(in), dimension(:) :: x, y, z
-    double precision s(size(x),size(y),size(z))
-    double precision, parameter :: pi = acos(-1D0)
+    real, intent(in), dimension(:) :: x, y, z
+    real s(size(x),size(y),size(z))
+    real, parameter :: pi = acos(-1E0)
     do concurrent(integer :: j=1:size(y), k=1:size(z)) default(none) shared(x,y,z,s)
       associate(eta => y(j)*cos(pi/4) + z(k)*sin(pi/4)) ! x-eta plane rotated around x axis pi/4 radians from x-y plane
         s(:,j,k) = (x**2 - eta**2)/2 
@@ -61,9 +61,9 @@ pure function rotated_stagnation_point_velocity(x,y,z) result(gradient)
     !! gradient(x,eta) = [ds/dx,              ds/dy,              ds/dz]
     !!                 = [    x, (ds/deta)(deta/dy), (ds/deta)(deta/dz)]
     !!                 = [    x,   -eta * cos(pi/4),   -eta * sin(pi/4)]
-    double precision, intent(in), dimension(:) :: x, y, z
-    double precision gradient(size(x),size(y),size(z),dimensionality)
-    double precision, parameter :: pi = acos(-1D0)
+    real, intent(in), dimension(:) :: x, y, z
+    real gradient(size(x),size(y),size(z),dimensionality)
+    real, parameter :: pi = acos(-1E0)
     do concurrent(integer :: i=1:size(x), j=1:size(y), k=1:size(z)) default(none) shared(gradient,x,y,z)
       associate(eta => y(j)*cos(pi/4) + z(k)*sin(pi/4)) ! x-eta plane rotated around x axis pi/4 radians from x-y plane
         gradient(i,j,k,:) = [x(i), -eta * cos(pi/4), -eta * sin(pi/4)]
@@ -82,7 +82,7 @@ function check_gradient() result(test_diagnosis)
     test_diagnosis = passing_test()
 
     do order = 2, 4, 2
-      associate(scalar_3D => scalar_3D_t(scalar_3D_initializer, order=order, cells=[20,20,20], x_min=[-10D0,-10D0,-10D0], x_max=[10D0,10D0,10D0]))
+      associate(scalar_3D => scalar_3D_t(scalar_3D_initializer, order=order, cells=[20,20,20], x_min=[-10E0,-10E0,-10E0], x_max=[10E0,10E0,10E0]))
         associate( &
            grad_scalar => .grad. scalar_3D &
           ,expected_gradient => vector_3D_t(expected_gradient_initializer, mold=scalar_3D) &
diff --git a/test/vector_2D_test_m.F90 b/test/vector_2D_test_m.F90
index 0ee6a60..2373f2e 100644
--- a/test/vector_2D_test_m.F90
+++ b/test/vector_2D_test_m.F90
@@ -35,8 +35,8 @@ module vector_2D_test_m
   end type
 
   integer, parameter :: space_dimension = 2
-  double precision, parameter :: tolerance = 1D-2
-  double precision, parameter :: u_const(*) = [1D0,2D0], v_const(*) = [3D0,4D0], u_dot_v_exact = dot_product(u_const, v_const)
+  real, parameter :: tolerance = 1E-2
+  real, parameter :: u_const(*) = [1E0,2E0], v_const(*) = [3E0,4E0], u_dot_v_exact = dot_product(u_const, v_const)
 
 contains
 
@@ -57,8 +57,8 @@ function results() result(test_results)
   end function
 
   pure function biquadratic(x,y) result(z)
-    double precision, intent(in) :: x(:), y(:)
-    double precision z(size(x),size(y),space_dimension)
+    real, intent(in) :: x(:), y(:)
+    real z(size(x),size(y),space_dimension)
     do concurrent(integer :: i=1:size(x), j=1:size(y)) default(none) shared(x,y,z)
       z(i,j,:) = [ &
          1 - 2*x(i) + 3*x(i)**2 - x(i)*y(j)/5 + 3*y(j)**2 - 2*y(j) &
@@ -68,16 +68,16 @@ pure function biquadratic(x,y) result(z)
   end function
 
   pure function biquadratic_divergence(x,y) result(divergence)
-    double precision, intent(in) :: x(:), y(:)
-    double precision divergence(size(x),size(y))
+    real, intent(in) :: x(:), y(:)
+    real divergence(size(x),size(y))
     do concurrent(integer :: i=1:size(x), j=1:size(y)) default(none) shared(divergence,x,y)
       divergence(i,j) = (-2 + 6*x(i) - y(j)/5) + (-2 + 6*y(j) - x(i)/5)
     end do
   end function
 
   pure function cubic(x,y) result(z)
-    double precision, intent(in) :: x(:), y(:)
-    double precision z(size(x),size(y),space_dimension)
+    real, intent(in) :: x(:), y(:)
+    real z(size(x),size(y),space_dimension)
     do concurrent(integer :: i=1:size(x), j=1:size(y)) default(none) shared(x,y,z)
       z(i,j,:) = [ &
          1 - 2*x(i) + 3*x(i)**3 - x(i)*y(j)/5 + 3*y(j)**3 - 2*y(j) &
@@ -87,8 +87,8 @@ pure function cubic(x,y) result(z)
   end function
 
   pure function cubic_divergence(x,y) result(divergence)
-    double precision, intent(in) :: x(:), y(:)
-    double precision divergence(size(x),size(y))
+    real, intent(in) :: x(:), y(:)
+    real divergence(size(x),size(y))
     do concurrent(integer :: i=1:size(x), j=1:size(y)) default(none) shared(divergence,x,y)
       divergence(i,j) = (-2 + 9*x(i)**2 - y(j)/5) + (-2 + 9*y(j)**2 - x(i)/5)
     end do
@@ -113,7 +113,7 @@ function check_divergence() result(test_diagnosis)
       case default
         error stop "check_divergence(vector_2D_test_m): unsupported order"
       end select
-      associate(vector_2D => vector_2D_t(vector_2D_initializer, order=order, cells=[40,30], x_min=[0D0,0D0], x_max=[2D0,1D0]))
+      associate(vector_2D => vector_2D_t(vector_2D_initializer, order=order, cells=[40,30], x_min=[0E0,0E0], x_max=[2E0,1E0]))
         associate(div_vector => .div. vector_2D)
           associate(expected_divergence => divergence_2D_t(expected_divergence_initializer, mold=vector_2D))
             test_diagnosis = test_diagnosis .also. &
@@ -126,16 +126,16 @@ function check_divergence() result(test_diagnosis)
   end function
 
   pure function u_field(x,y) result(u)
-    double precision, intent(in) :: x(:), y(:)
-    double precision u(size(x),size(y),space_dimension)
+    real, intent(in) :: x(:), y(:)
+    real u(size(x),size(y),space_dimension)
     do concurrent(integer :: i=1:size(x), j=1:size(y))
       u(i,j,:) = u_const
     end do
   end function
 
   pure function v_field(x,y) result(v)
-    double precision, intent(in) :: x(:), y(:)
-    double precision v(size(x),size(y),space_dimension)
+    real, intent(in) :: x(:), y(:)
+    real v(size(x),size(y),space_dimension)
     do concurrent(integer :: i=1:size(x), j=1:size(y))
       v(i,j,:) = v_const
     end do
@@ -153,35 +153,35 @@ function check_dot_product() result(test_diagnosis)
 
     do order = 2, 4, 2
       associate( &
-         u => vector_2D_t(u_init, cells=[10,10], x_min=[0D0,0D0], x_max=[5D0,5D0], order=order) &
-        ,v => vector_2D_t(v_init, cells=[10,10], x_min=[0D0,0D0], x_max=[5D0,5D0], order=order) &
+         u => vector_2D_t(u_init, cells=[10,10], x_min=[0E0,0E0], x_max=[5E0,5E0], order=order) &
+        ,v => vector_2D_t(v_init, cells=[10,10], x_min=[0E0,0E0], x_max=[5E0,5E0], order=order) &
       )
         associate(u_dot_v => u .dot. v)
-          test_diagnosis = test_diagnosis .also. (.all. (u_dot_v%values() .approximates. u_dot_v_exact .within. 1D-6))
+          test_diagnosis = test_diagnosis .also. (.all. (u_dot_v%values() .approximates. u_dot_v_exact .within. 1E-5))
         end associate
       end associate
     end do
   end function
 
   pure function scalar_field(x,y) result(s)
-    double precision, intent(in) :: x(:), y(:)
-    double precision s(size(x),size(y))
+    real, intent(in) :: x(:), y(:)
+    real s(size(x),size(y))
     do concurrent(integer :: i=1:size(x), j=1:size(y))
       s(i,j) = x(i)
     end do
   end function
 
   pure function vector_field(x,y) result(v)
-    double precision, intent(in) :: x(:), y(:)
-    double precision v(size(x),size(y),space_dimension)
+    real, intent(in) :: x(:), y(:)
+    real v(size(x),size(y),space_dimension)
     do concurrent(integer :: i=1:size(x), j=1:size(y))
       v(i,j,:) = y(j)
     end do
   end function
 
   pure function scalar_vector_product(x,y) result(u)
-    double precision, intent(in) :: x(:), y(:)
-    double precision u(size(x),size(y),space_dimension)
+    real, intent(in) :: x(:), y(:)
+    real u(size(x),size(y),space_dimension)
     do concurrent(integer :: i=1:size(x), j=1:size(y))
       u(i,j,:) = x(i)*y(j)
     end do
@@ -201,13 +201,13 @@ function check_vector_scalar_product() result(test_diagnosis)
 
     do order = 2, 4, 2
       associate( &
-         s => scalar_2D_t(s_init, cells=[10,10], x_min=[0D0,0D0], x_max=[5D0,5D0], order=order) &
-        ,v => vector_2D_t(v_init, cells=[10,10], x_min=[0D0,0D0], x_max=[5D0,5D0], order=order) &
-        ,vs_expected => vector_2D_t(vs_init, cells=[10,10], x_min=[0D0,0D0], x_max=[5D0,5D0], order=order) &
+         s => scalar_2D_t(s_init, cells=[10,10], x_min=[0E0,0E0], x_max=[5E0,5E0], order=order) &
+        ,v => vector_2D_t(v_init, cells=[10,10], x_min=[0E0,0E0], x_max=[5E0,5E0], order=order) &
+        ,vs_expected => vector_2D_t(vs_init, cells=[10,10], x_min=[0E0,0E0], x_max=[5E0,5E0], order=order) &
       )
         associate(vs => v * s)
-          test_diagnosis = test_diagnosis .also. (.all. (vs%values(x_dir) .approximates. vs_expected%values(x_dir) .within. 1D-6))
-          test_diagnosis = test_diagnosis .also. (.all. (vs%values(y_dir) .approximates. vs_expected%values(y_dir) .within. 1D-6))
+          test_diagnosis = test_diagnosis .also. (.all. (vs%values(x_dir) .approximates. vs_expected%values(x_dir) .within. 1E-5))
+          test_diagnosis = test_diagnosis .also. (.all. (vs%values(y_dir) .approximates. vs_expected%values(y_dir) .within. 1E-5))
         end associate
       end associate
     end do
diff --git a/test/vector_3D_test_m.F90 b/test/vector_3D_test_m.F90
index 5e39137..11ac916 100644
--- a/test/vector_3D_test_m.F90
+++ b/test/vector_3D_test_m.F90
@@ -38,7 +38,7 @@ module vector_3D_test_m
     procedure, nopass :: results
   end type
 
-  double precision, parameter :: tolerance = 1D-12
+  real, parameter :: tolerance = 1E-4
 
 contains
 
@@ -60,9 +60,9 @@ function results() result(test_results)
 
   pure function rotated_stagnation_point_potential(x,y,z) result(s)
     !! Define a stagnation-point scalar potential in a plane tilted 45 deg from the x-y plane
-    double precision, intent(in), dimension(:) :: x, y, z
-    double precision s(size(x),size(y),size(z))
-    double precision, parameter :: pi = acos(-1D0)
+    real, intent(in), dimension(:) :: x, y, z
+    real s(size(x),size(y),size(z))
+    real, parameter :: pi = acos(-1E0)
     do concurrent(integer :: j=1:size(y), k=1:size(z)) default(none) shared(x,y,z,s)
       associate(eta => y(j)*cos(pi/4) + z(k)*sin(pi/4)) ! x-eta plane rotated around x axis pi/4 radians from x-y plane
         s(:,j,k) = (x**2 - eta**2)/2
@@ -75,10 +75,10 @@ pure function rotated_stagnation_point_velocity(x,y,z) result(gradient)
     !! gradient(x,eta) = [ds/dx,              ds/dy,              ds/dz]
     !!                 = [    x, (ds/deta)(deta/dy), (ds/deta)(deta/dz)]
     !!                 = [    x,  -eta * cos(theta),  -eta * sin(theta)]
-    double precision, intent(in), dimension(:) :: x, y, z
-    double precision, parameter :: pi = acos(-1D0), theta = pi/4
+    real, intent(in), dimension(:) :: x, y, z
+    real, parameter :: pi = acos(-1E0), theta = pi/4
     integer, parameter :: dimensionality = 3
-    double precision gradient(size(x),size(y),size(z),dimensionality)
+    real gradient(size(x),size(y),size(z),dimensionality)
     do concurrent(integer :: i=1:size(x), j=1:size(y), k=1:size(z)) default(none) shared(gradient,x,y,z)
       associate(eta => y(j)*cos(theta) + z(k)*sin(theta)) ! x-eta plane rotated around x axis theta radians from x-y plane
         gradient(i,j,k,:) = [x(i), -eta * cos(theta), -eta * sin(theta)]
@@ -87,10 +87,10 @@ pure function rotated_stagnation_point_velocity(x,y,z) result(gradient)
   end function
 
   pure function velocity_squared(x,y,z) result(v_sq)
-    double precision, intent(in), dimension(:) :: x, y, z
-    double precision, parameter :: pi = acos(-1D0), theta = pi/4
+    real, intent(in), dimension(:) :: x, y, z
+    real, parameter :: pi = acos(-1E0), theta = pi/4
     integer, parameter :: dimensionality = 3
-    double precision v_sq(size(x),size(y),size(z))
+    real v_sq(size(x),size(y),size(z))
     do concurrent(integer :: i=1:size(x), j=1:size(y), k=1:size(z)) default(none) shared(v_sq, x, y, z)
       associate(eta => y(j)*cos(theta) + z(k)*sin(theta)) ! x-eta plane rotated around x axis theta radians from x-y plane
         associate(v => [x(i), -eta * cos(theta), -eta * sin(theta)])
@@ -111,11 +111,11 @@ function check_divergence() result(test_diagnosis)
     !expected_divergence_initializer => cubic_divergence
 
     do order = 2, 4, 2
-      associate(vector_3D => vector_3D_t(vector_3D_initializer, order, cells=[20,20,20], x_min=[-100D0,-10D0,-10D0], x_max=[10D0,10D0,10D0]))
+      associate(vector_3D => vector_3D_t(vector_3D_initializer, order, cells=[20,20,20], x_min=[-100E0,-10E0,-10E0], x_max=[10E0,10E0,10E0]))
         associate(div_vector => .div. vector_3D)
           !associate(expected_divergence => divergence_3D_t(expected_divergence_initializer, mold=vector_3D))
             test_diagnosis = test_diagnosis .also. &
-              (.all. (div_vector%values() .approximates. 0D0 .within. tolerance)) &
+              (.all. (div_vector%values() .approximates. 0E0 .within. tolerance)) &
               // string_t(" for order ") // string_t(order)
           !end associate
         end associate
@@ -136,8 +136,8 @@ function check_dot_product() result(test_diagnosis)
 
     do order = 2, 4, 2
       associate( &
-         v => vector_3D_t(v_init, order, cells=[20,20,20], x_min=[-10D0, -10D0, -10D0], x_max=[10D0, 10D0, 10D0]) &
-        ,v_squared => scalar_3D_t(v_sq_init, order, cells=[20,20,20], x_min=[-10D0, -10D0, -10D0], x_max=[10D0, 10D0, 10D0]) &
+         v => vector_3D_t(v_init, order, cells=[20,20,20], x_min=[-10E0, -10E0, -10E0], x_max=[10E0, 10E0, 10E0]) &
+        ,v_squared => scalar_3D_t(v_sq_init, order, cells=[20,20,20], x_min=[-10E0, -10E0, -10E0], x_max=[10E0, 10E0, 10E0]) &
       )
         associate(v_dot_v => v .dot. v)
           test_diagnosis = test_diagnosis .also. (.all. (v_dot_v%values() .approximates. v_squared%values() .within. tolerance))
@@ -158,7 +158,7 @@ function check_vector_scalar_product() result(test_diagnosis)
     s_init => rotated_stagnation_point_potential
 
     do order = 2, 4, 2
-      associate(s => scalar_3D_t(s_init, order, cells=[10,10,10], x_min=[-10D0, -10D0, -10D0], x_max=[10D0, 10D0, 10D0]))
+      associate(s => scalar_3D_t(s_init, order, cells=[10,10,10], x_min=[-10E0, -10E0, -10E0], x_max=[10E0, 10E0, 10E0]))
         associate(v => vector_3D_t(v_init, mold = s))
           associate(vs => v * s)
             test_diagnosis = test_diagnosis .also. (.all. (vs%values(x_dir) .approximates. v%values(x_dir)*s%to_faces(x_dir) .within. tolerance))