diff --git a/fortitude.toml b/fortitude.toml
new file mode 100644
index 00000000..a2e923b6
--- /dev/null
+++ b/fortitude.toml
@@ -0,0 +1,12 @@
+[check]
+ignore = [
+  "S001", # line-too-long
+  "S051", # deprecated-relational-operator (fix available)
+  "S101", # trailing-whitespace
+  "S102", # incorrect-space-before-comment, we should enable this, autofix enable
+  "M011", # use-all (should use "only:")
+  "M001", # procedure-not-in-module (need a per-file ignore for init.F90)
+  "T042", # assumed-size-character-intent TODO: We should fix all instances!
+]
+line-length = 132
+show-fixes = true
diff --git a/src/analysis.F90 b/src/analysis.F90
index 20400f0a..f95105e9 100644
--- a/src/analysis.F90
+++ b/src/analysis.F90
@@ -28,7 +28,6 @@ subroutine analysis(x, y, z, vx, vy, vz, fxc, fyc, fzc, eclas)
       use mod_analyze_geometry
       use mod_io
       use mod_system, only: am
-      implicit none
       real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(in) :: fxc(:, :), fyc(:, :), fzc(:, :)
       real(DP), intent(in) :: vx(:, :), vy(:, :), vz(:, :)
@@ -102,7 +101,6 @@ subroutine trajout(x, y, z, time_step)
       use mod_general, only: nwalk, natom, sim_time
       use mod_mpi, only: get_mpi_rank
       use mod_system, only: names
-      implicit none
       real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
       integer, intent(in) :: time_step
       integer :: iat, iw
diff --git a/src/cmdline.F90 b/src/cmdline.F90
index 3b1ac40b..30473611 100644
--- a/src/cmdline.F90
+++ b/src/cmdline.F90
@@ -1,12 +1,12 @@
 module mod_cmdline
    use mod_interfaces, only: print_compile_info
+   implicit none
    private
    public :: get_cmdline
 
 contains
 
    subroutine print_help()
-      implicit none
       integer, dimension(8) :: time_data
 
       call date_and_time(values=time_data)
diff --git a/src/ekin.F90 b/src/ekin.F90
index 7beef370..c74e04c4 100644
--- a/src/ekin.F90
+++ b/src/ekin.F90
@@ -15,7 +15,6 @@ subroutine temperature(px, py, pz, amt, eclas)
       use mod_nhc, only: inose, get_nhcham
       use mod_gle, only: get_langham
       use mod_shake, only: nshake
-      implicit none
       integer :: iw, iat
       real(DP), intent(in) :: px(:, :), py(:, :), pz(:, :)
       real(DP), intent(in) :: amt(:, :)
@@ -68,7 +67,6 @@ end subroutine temperature
    real(DP) function ekin_v(vx, vy, vz)
       use mod_general, only: nwalk, natom
       use mod_system, only: am
-      implicit none
       real(DP), intent(in) :: vx(:, :), vy(:, :), vz(:, :)
       real(DP) :: temp1, ekin_mom
       integer :: iw, iat
diff --git a/src/en_restraint.F90 b/src/en_restraint.F90
index 14be4820..dde43bee 100644
--- a/src/en_restraint.F90
+++ b/src/en_restraint.F90
@@ -109,7 +109,7 @@ subroutine energy_restraint(x, y, z, px, py, pz, eclas)
 
             !Iterative procedure
             convercrit = 100
-            do while (convercrit > 0.00001)
+            do while (convercrit > 0.00001D0)
                ! deltaE(t+dt) prediction - improves the accuracy
                deltaEnext = 0
                do iat = 1, natom
@@ -166,7 +166,7 @@ subroutine energy_restraint(x, y, z, px, py, pz, eclas)
             eclas = eclas ! + quadratic_restraint_energy
 
             !Output to en_restraint.dat
-            write (UERMD, '(I8,F16.8,E20.10,E20.10,F16.8)') it, excE * AUTOEV, deltaE, 0.0, 0.0
+            write (UERMD, '(I8,F16.8,E20.10,E20.10,F16.8)') it, excE * AUTOEV, deltaE, 0.0D0, 0.0D0
 
          end if
 
diff --git a/src/error.F90 b/src/error.F90
index 6d0ee452..b232bffb 100644
--- a/src/error.F90
+++ b/src/error.F90
@@ -28,6 +28,7 @@ module mod_error
 
    abstract interface
       subroutine error(filename, line_number, message)
+         implicit none
          character(len=*), intent(in) :: filename
          integer, intent(in) :: line_number
          character(len=*), intent(in) :: message
diff --git a/src/fftw_interface.F90 b/src/fftw_interface.F90
index ab90fe17..b782f5ec 100644
--- a/src/fftw_interface.F90
+++ b/src/fftw_interface.F90
@@ -38,7 +38,7 @@ subroutine fftw_normalmodes_init(nwalk)
 
       deallocate (x_tmp)
       deallocate (cx_tmp)
-   end subroutine
+   end subroutine fftw_normalmodes_init
 
    ! Simple wrapper functions around the FFTW interface.
    subroutine dft_normalmode2cart(nm, cart)
@@ -72,16 +72,16 @@ subroutine fftw_normalmodes_init(nwalk)
    end subroutine fftw_normalmodes_init
 
    subroutine dft_normalmode2cart(nm, cart)
-      complex(C_DOUBLE_COMPLEX), dimension(:) :: nm
-      real(C_DOUBLE), dimension(:) :: cart
+      complex(C_DOUBLE_COMPLEX), dimension(:), intent(inout) :: nm
+      real(C_DOUBLE), dimension(:), intent(inout) :: cart
       cart = 0.0D0
       nm = (0.0D0, 0.0D0)
       call not_compiled_with('FFTW library')
    end subroutine dft_normalmode2cart
 
    subroutine dft_cart2normalmode(cart, nm)
-      complex(C_DOUBLE_COMPLEX), dimension(:) :: nm
-      real(C_DOUBLE), dimension(:) :: cart
+      complex(C_DOUBLE_COMPLEX), dimension(:), intent(inout) :: nm
+      real(C_DOUBLE), dimension(:), intent(inout) :: cart
       cart = 0.0D0
       nm = (0.0D0, 0.0D0)
       call not_compiled_with('FFTW library')
diff --git a/src/force_abin.F90 b/src/force_abin.F90
index 58c5cdcb..82ac798f 100644
--- a/src/force_abin.F90
+++ b/src/force_abin.F90
@@ -165,7 +165,7 @@ real(DP) function read_energy(engrad_unit, abort) result(energy)
 !$OMP FLUSH(abort)
          return
       end if
-   end function
+   end function read_energy
 
    subroutine read_forces(fx, fy, fz, num_atom, iw, engrad_unit, abort)
       use mod_files, only: stderr
@@ -209,7 +209,6 @@ subroutine force_abin(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
       use mod_lz, only: nstate_lz, tocalc_lz, en_array_lz, istate_lz, write_lz_data
       use mod_qmmm, only: natqm
       use mod_utils, only: toupper, append_rank
-      implicit none
       real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(out) :: fx(:, :), fy(:, :), fz(:, :)
       real(DP), intent(out) :: eclas
@@ -325,7 +324,6 @@ subroutine oniom(x, y, z, fx, fy, fz, eclas, iw, abort)
       use mod_qmmm, only: natqm
       use mod_sh_integ, only: nstate
       use mod_sh, only: en_array, istate
-      implicit none
       real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :)
       real(DP), intent(inout) :: eclas
diff --git a/src/force_bound.F90 b/src/force_bound.F90
index d9911b90..ab4aa047 100644
--- a/src/force_bound.F90
+++ b/src/force_bound.F90
@@ -61,7 +61,7 @@ subroutine sbc_init(x, y, z)
          write (*, *) 'Calculating cluster radius from given densty.'
          rho = rho * fact !conversion from g/L to atomic units
          rb_sbc = mass_total / rho * 3 / 4 / PI
-         rb_sbc = rb_sbc**(1 / 3.)
+         rb_sbc = rb_sbc**(1 / 3.0D0)
       end if
 
       if (rmax > rb_sbc) then
@@ -73,7 +73,7 @@ subroutine sbc_init(x, y, z)
 
       write (*, *) 'rb_sbc[A]=', rb_sbc / ang
 
-   end subroutine
+   end subroutine sbc_init
 
    subroutine force_sbc(x, y, z, fx, fy, fz, walkmax)
       use mod_const, only: ANG
@@ -118,7 +118,7 @@ subroutine force_sbc(x, y, z, fx, fy, fz, walkmax)
       end if
 
       return
-   end subroutine
+   end subroutine force_sbc
 
    ! TODO: This could be a general purpose routine,
    ! move to utils
diff --git a/src/force_h2o.F90 b/src/force_h2o.F90
index e8ddba68..733332e8 100644
--- a/src/force_h2o.F90
+++ b/src/force_h2o.F90
@@ -129,6 +129,13 @@ subroutine force_h2o_cvrqd(x, y, z, fx, fy, fz, Eclas, natom, nbeads)
       ! TODO: Given the small difference between the Schwenke potential,
       ! we might not need to implement numerical forces here.
       ! call numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads)
+      ! Just to squash compiler warnings
+      if (.false.) then
+         print*,natom
+         fx = 0.0D0
+         fy = 0.0D0
+         fz = 0.0D0
+      end if
 
    end subroutine force_h2o_cvrqd
 
@@ -159,7 +166,7 @@ subroutine numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads)
       real(DP) :: Epot_delta(1)
 
       real(DP) :: Eclas_orig
-      real(DP) :: delta = 5.0E-5_DP
+      real(DP), parameter :: DELTA = 5.0E-5_DP
       integer :: i, j, k
 
       ! Calculate forces numerically using central differences
@@ -180,11 +187,11 @@ subroutine numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads)
                ! Move the atom forwards
                select case (k)
                case (1)
-                  x_new_forward(i, j) = x_new_forward(i, j) + delta
+                  x_new_forward(i, j) = x_new_forward(i, j) + DELTA
                case (2)
-                  y_new_forward(i, j) = y_new_forward(i, j) + delta
+                  y_new_forward(i, j) = y_new_forward(i, j) + DELTA
                case (3)
-                  z_new_forward(i, j) = z_new_forward(i, j) + delta
+                  z_new_forward(i, j) = z_new_forward(i, j) + DELTA
                end select
 
                ! Calculate the energy for the forward perturbed geometry
@@ -197,11 +204,11 @@ subroutine numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads)
                ! Calculate the numerical force
                select case (k)
                case (1)
-                  fx(i, j) = -(Epot_delta(1) - Eclas_orig) / delta
+                  fx(i, j) = -(Epot_delta(1) - Eclas_orig) / DELTA
                case (2)
-                  fy(i, j) = -(Epot_delta(1) - Eclas_orig) / delta
+                  fy(i, j) = -(Epot_delta(1) - Eclas_orig) / DELTA
                case (3)
-                  fz(i, j) = -(Epot_delta(1) - Eclas_orig) / delta
+                  fz(i, j) = -(Epot_delta(1) - Eclas_orig) / DELTA
                end select
             end do
          end do
diff --git a/src/force_spline.F90 b/src/force_spline.F90
index 88684f08..394edc01 100644
--- a/src/force_spline.F90
+++ b/src/force_spline.F90
@@ -87,7 +87,7 @@ subroutine finalize_spline()
       if (allocated(second_derivatives)) then
          deallocate (second_derivatives)
       end if
-   end subroutine
+   end subroutine finalize_spline
 
    subroutine read_grid(fname, x_grid, y_grid, grid_size)
       use mod_files, only: stdout
@@ -122,7 +122,7 @@ subroutine read_grid(fname, x_grid, y_grid, grid_size)
       end do
       close (u)
 
-   end subroutine
+   end subroutine read_grid
 
    subroutine validate_grid(x_grid, ngrid)
       real(DP), dimension(ngrid), intent(in) :: x_grid
@@ -141,7 +141,7 @@ subroutine validate_grid(x_grid, ngrid)
             return
          end if
       end do
-   end subroutine
+   end subroutine validate_grid
 
    subroutine print_splined_potential(fname, x_grid, grid_size)
       character(len=*), intent(in) :: fname
@@ -158,6 +158,6 @@ subroutine print_splined_potential(fname, x_grid, grid_size)
          x = x + dx
       end do
       close (u)
-   end subroutine
+   end subroutine print_splined_potential
 
-end module
+end module mod_splined_grid
diff --git a/src/force_tcpb.F90 b/src/force_tcpb.F90
index 8acf5328..b0603b82 100644
--- a/src/force_tcpb.F90
+++ b/src/force_tcpb.F90
@@ -33,7 +33,7 @@ module mod_force_tcpb
       character(len=1024) :: input_file = ''
       ! This make TC reuse WF from previous step.
       integer :: globaltreatment = 0
-   end type
+   end type tcpb_params
    type(tcpb_params) :: tcpb
    character(len=5), allocatable :: qmattypes(:)
    save
diff --git a/src/forces.F90 b/src/forces.F90
index 5af786cb..a04d3f9e 100644
--- a/src/forces.F90
+++ b/src/forces.F90
@@ -260,9 +260,9 @@ subroutine force_quantum(fx, fy, fz, x, y, z, amg, quantum_energy)
             fz(j, i) = (z(j, i) - z(j, kplus))
             fz(j, i) = fz(j, i) + (z(j, i) - z(j, kminus))
             fz(j, i) = -fz(j, i) * ak(j, i)
-            equant = equant + 0.5 * ak(j, i) * (x(j, i) - x(j, kplus))**2
-            equant = equant + 0.5 * ak(j, i) * (y(j, i) - y(j, kplus))**2
-            equant = equant + 0.5 * ak(j, i) * (z(j, i) - z(j, kplus))**2
+            equant = equant + 0.5D0 * ak(j, i) * (x(j, i) - x(j, kplus))**2
+            equant = equant + 0.5D0 * ak(j, i) * (y(j, i) - y(j, kplus))**2
+            equant = equant + 0.5D0 * ak(j, i) * (z(j, i) - z(j, kplus))**2
          end do
       end do
    end if
diff --git a/src/fortran_interfaces.F90 b/src/fortran_interfaces.F90
index 5648a7e2..182937a3 100644
--- a/src/fortran_interfaces.F90
+++ b/src/fortran_interfaces.F90
@@ -6,18 +6,22 @@
 module mod_interfaces
    use, intrinsic :: iso_c_binding, only: C_INT, C_INT32_T
    use mod_const, only: DP
+   implicit none
    public
    interface
 
       subroutine print_compile_info()
+         implicit none
       end subroutine print_compile_info
 
       subroutine finish(error_code)
+         implicit none
          integer, intent(in) :: error_code
       end subroutine finish
 
       subroutine force_clas(fx, fy, fz, x, y, z, eclas, chpot)
          import :: DP
+         implicit none
          real(DP), intent(inout) :: x(:, :), y(:, :), z(:, :)
          real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :)
          real(DP), intent(out) :: eclas
@@ -26,6 +30,7 @@ end subroutine force_clas
 
       subroutine force_quantum(fx, fy, fz, x, y, z, amg, energy)
          import :: DP
+         implicit none
          real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
          real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :)
          real(DP), intent(in) :: amg(:, :)
@@ -34,6 +39,7 @@ end subroutine force_quantum
 
       subroutine force_wrapper(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
          import :: DP
+         implicit none
          real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
          real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :)
          real(DP), intent(out) :: eclas
@@ -42,6 +48,7 @@ subroutine force_wrapper(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
       end subroutine force_wrapper
 
       subroutine omp_set_num_threads(nthreads)
+         implicit none
          integer, intent(in) :: nthreads
       end subroutine omp_set_num_threads
 
@@ -49,6 +56,7 @@ end subroutine omp_set_num_threads
       ! of how we pass the 2D arrays...
       !subroutine force_water(x, y, z, fx, fy, fz, eclas, natom, walkmax, watpot) !bind(C)
       !   use, intrinsic :: iso_c_binding
+      !   implicit none
       !   real(C_DOUBLE), intent(in) :: x(:, :), y(:, :), z(:, :)
       !   real(C_DOUBLE), intent(inout) :: fx(:, :), fy(:, :), fz(:, :)
       !   real(C_DOUBLE), intent(inout) :: eclas
@@ -59,7 +67,8 @@ end subroutine omp_set_num_threads
       ! int usleep(useconds_t useconds)
       function usleep(useconds) bind(c, name='usleep')
          import :: C_INT, C_INT32_T
-         integer(kind=C_INT32_T), value :: useconds
+         implicit none
+         integer(kind=C_INT32_T), intent(in), value :: useconds
          integer(kind=C_INT) :: usleep
       end function usleep
 
@@ -67,6 +76,7 @@ end function usleep
       ! using Schwenke potential, see h2o_schwenke.f
       subroutine h2o_pot_schwenke(rij, v, n)
          import :: DP
+         implicit none
          integer, intent(in) :: n
          real(DP), intent(in) :: rij(n, 3)
          real(DP), intent(out) :: v(n)
@@ -74,6 +84,7 @@ end subroutine h2o_pot_schwenke
 
       subroutine h2o_pot_cvrqd(V, rOH1, rOH2, aHOH, mH, mO)
          import :: DP
+         implicit none
          real(DP), intent(out) :: V
          real(DP), intent(in) :: rOH1, rOH2, aHOH, mH, mO
       end subroutine h2o_pot_cvrqd
diff --git a/src/gle.F90 b/src/gle.F90
index 5173456b..e8a4bb83 100644
--- a/src/gle.F90
+++ b/src/gle.F90
@@ -48,7 +48,7 @@ module mod_gle_private
    ! Used both by white-noise and GLE thermostats.
    real(DP) function get_langham()
       get_langham = langham
-   end function
+   end function get_langham
 
    ! Initialize white-noise PILE thermostat,
    ! which can be used both for classical MD and PIMD.
@@ -93,7 +93,7 @@ subroutine pile_init(dt, tau0)
          c1(iw) = exp(-dt * gam)
          c2(iw) = dsqrt(1 - c1(iw)**2) * dsqrt(temp * nwalk)
       end do
-   end subroutine
+   end subroutine pile_init
 
    ! White-noise propagator. time-step has been set in pile_init
    subroutine pile_step(px, py, pz, m)
@@ -124,7 +124,7 @@ subroutine pile_step(px, py, pz, m)
       end do
 
       langham = langham - ekin_p(px, py, pz, m, natom, nwalk)
-   end subroutine
+   end subroutine pile_step
 
    subroutine finalize_pile()
       if (allocated(c1)) deallocate (c1, c2)
@@ -222,7 +222,6 @@ subroutine gle_init(dt)
       use mod_mpi, only: get_mpi_rank
       use mod_general, only: natom, nwalk, ipimd, inormalmodes
       use mod_nhc, only: temp, inose
-      implicit none
       real(DP), intent(in) :: dt
       real(DP), allocatable :: gA(:, :), gC(:, :)
       character(len=100) :: glea, glec, glea_centroid, glec_centroid
@@ -351,12 +350,12 @@ subroutine pile_restin(u)
       ! Restart file unit
       integer, intent(in) :: u
       read (u, *) langham
-   end subroutine
+   end subroutine pile_restin
 
    subroutine pile_restout(u)
       integer, intent(in) :: u
       write (u, '(ES24.16E3)') langham
-   end subroutine
+   end subroutine pile_restout
 
    subroutine gle_restin(u)
       use mod_general, only: natom, nwalk
@@ -373,7 +372,7 @@ subroutine gle_restin(u)
          end do
       end do
       read (u, *) langham
-   end subroutine
+   end subroutine gle_restin
 
    subroutine gle_restout(u)
       use mod_general, only: natom, nwalk
@@ -389,7 +388,7 @@ subroutine gle_restout(u)
          end do
       end do
       write (u, *) langham
-   end subroutine
+   end subroutine gle_restout
 
    subroutine finalize_gle()
       if (allocated(gS)) then
@@ -549,7 +548,7 @@ subroutine matrix_exp(M, n, j, k, EM)
       real(DP), intent(out) :: EM(n, n)
 
       real(DP) :: tc(j + 1), SM(n, n)
-      integer p, i
+      integer :: p, i
       tc(1) = 1
       do i = 1, j
          tc(i + 1) = tc(i) / dble(i)
@@ -586,7 +585,7 @@ subroutine cholesky(SST, S, n)
       real(DP), intent(in) :: SST(n, n)
       real(DP), intent(out) :: S(n, n)
       real(DP) :: L(n, n), D(n, n)
-      integer i, j, k
+      integer :: i, j, k
       S = 0.D0
       L = 0.D0
       D = 0.D0
diff --git a/src/init.F90 b/src/init.F90
index 92b64426..0804fab8 100644
--- a/src/init.F90
+++ b/src/init.F90
@@ -67,13 +67,13 @@ subroutine init(dt)
       use mod_terampi
       use mod_terampi_sh
       use mod_mdstep, only: initialize_integrator, nabin, nstep_ref
-      implicit none
       real(DP), intent(out) :: dt
       ! Input parameters for analytical potentials
-      real(DP) :: lambda_dw = -1.0D0, D0_dw = -1.0D0, k_dw = -1.0D0, r0_dw = -1.D0
-      real(DP) :: r0_morse = -1, d0_morse = -1, k_morse = -1
-      real(DP) :: k = -1, r0 = -1
-      real(DP) :: kx = -1, ky = -1, kz = -1
+      ! TODO: Initialize these variable in the code not here
+      real(DP), save :: lambda_dw = -1.0D0, D0_dw = -1.0D0, k_dw = -1.0D0, r0_dw = -1.D0
+      real(DP), save :: r0_morse = -1, d0_morse = -1, k_morse = -1
+      real(DP), save :: k = -1, r0 = -1
+      real(DP), save :: kx = -1, ky = -1, kz = -1
       ! Lennard-Jones parameteres and Coulomb charges (pot=_mm_)
       ! All input parameters are expected to be in atomic units,
       ! except LJ_rmin which should be in angstroms.
@@ -84,14 +84,14 @@ subroutine init(dt)
       ! L-J parameters
       real(DP), allocatable :: LJ_rmin(:), LJ_eps(:)
       ! Initial temperature (read from namelist nhcopt)
-      real(DP) :: temp0 = -1
+      real(DP), save :: temp0 = -1
       ! User-defined masses in relative atomic units
       real(DP), allocatable :: masses(:)
       integer :: iw, iat, natom_xyz, iost
       integer :: shiftdihed
       ! Random number seed
       ! Negative value means we get the seed from /dev/urandom
-      integer :: irandom = -1
+      integer, save :: irandom = -1
       ! Number of OpenMP processes, read from ABIN input
       ! WARNING: We do NOT use OMP_NUM_THREADS environment variable!
       integer :: nproc
@@ -433,11 +433,11 @@ subroutine init(dt)
       if (inose == 1) then
          call nhc_init()
       else if (inose == 2) then
-         call gle_init(dt * 0.5 / nabin / nstep_ref) !nabin is set to 1 unless ipimd=1
+         call gle_init(dt * 0.5D0 / nabin / nstep_ref) !nabin is set to 1 unless ipimd=1
       else if (inose == 3) then
-         call pile_init(dt * 0.5, tau0_langevin)
+         call pile_init(dt * 0.5D0, tau0_langevin)
       else if (inose == 4) then
-         call gle_init(dt * 0.5 / nstep_ref)
+         call gle_init(dt * 0.5D0 / nstep_ref)
       else if (inose == 0) then
          write (stdout, '(A)') 'No thermostat. NVE ensemble.'
       else
diff --git a/src/landau_zener.F90 b/src/landau_zener.F90
index c8f4f4cc..cd2dddab 100644
--- a/src/landau_zener.F90
+++ b/src/landau_zener.F90
@@ -213,7 +213,7 @@ subroutine lz_hop(x, y, z, vx, vy, vz, fxc, fyc, fzc, amt, dt, eclas, chpot)
       integer :: S_to_T !itest, iost
       integer :: ist ! =istate_lz
       real(DP) :: Ekin, Epot, Epot2, dE, hop_rdnum, hop_rdnum2, vel_rescale, stepfs, molveloc, grad_diff
-      real(DP) :: one = 1.0D0
+      real(DP), parameter :: one = 1.0D0
       character(len=100) :: formt, fmt_in, fmt_out
 
       !TODO: energy drift
@@ -268,7 +268,7 @@ subroutine lz_hop(x, y, z, vx, vy, vz, fxc, fyc, fzc, amt, dt, eclas, chpot)
             ! Adding check for false 3P-minima
             ! only for significant probabilities, for tiny probabilities this does not make sense
             ! e.g. for parallel states (the whole LZ does not make sense for this case)
-            if (prob(ist1) > 0.01) then
+            if (prob(ist1) > 0.01D0) then
                ! Calculating backward second derivative formula. We are not using the newest energy
                ! but the previous three. Discontinuity would not affect this formula.
                second_der_back = ((en_diff(2) - 2 * en_diff(3) + en_diff(4)) / dt**2)
@@ -276,14 +276,14 @@ subroutine lz_hop(x, y, z, vx, vy, vz, fxc, fyc, fzc, amt, dt, eclas, chpot)
                ! If the change is too large, we either almost hit the CI or we have discontinuity.
                der_check = abs((second_der - second_der_back) / second_der)
                ! If they differ by more then 130%, we have aa unphysical change of curvature --> certain discontinuity.
-               if (der_check > 1.3) then
+               if (der_check > 1.3D0) then
                   write (stdout, *) "ERROR: Change of curvature --> discontinuity in PES!"
                   write (stdout, *) "Probability set to 0!"
                   prob(ist1) = 0.0D0
                   ! 30% threshold was set empirically and should capture most discontinuities
                   ! yet it can also be a conical intersection. Thus, we just issue an warning
                   ! and let the user to evaluate on his own.
-               else if (der_check > 0.3) then
+               else if (der_check > 0.3D0) then
                   write (stdout, *) "WARNING: Possible discontinuity in PES! Check PES.dat!"
                end if
             end if
@@ -605,9 +605,11 @@ subroutine lz_getgraddiff(grad_diff, ist, ist1, x, y, z, vx, vy, vz, fxc, fyc, f
       real(DP), intent(inout) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(inout) :: vx(:, :), vy(:, :), vz(:, :)
       real(DP), intent(inout) :: fxc(:, :), fyc(:, :), fzc(:, :)
-      real(DP) :: eclas = 0.0D0
+      real(DP) :: eclas
       integer :: i, iat
 
+      eclas = 0.0D0
+
       !Gradient to compute
       do i = 1, nstate_lz
          if (i == ist1) then
diff --git a/src/mdstep.F90 b/src/mdstep.F90
index 33fe85c0..8e7f7760 100644
--- a/src/mdstep.F90
+++ b/src/mdstep.F90
@@ -20,6 +20,7 @@ module mod_mdstep
    abstract interface
       subroutine integrator(x, y, z, px, py, pz, amt, dt, E_pot, fx, fy, fz)
          import :: DP
+         implicit none
          real(DP), dimension(:, :), intent(inout) :: x, y, z, px, py, pz
          real(DP), dimension(:, :), intent(in) :: amt
          real(DP), intent(in) :: dt
@@ -44,7 +45,7 @@ subroutine initialize_integrator(dt, ipimd, inormalmodes, nshake, pot, pot_ref)
       use mod_files, only: stdout
       real(DP), intent(in) :: dt
       integer, intent(in) :: ipimd, inormalmodes, nshake
-      character(len=*) :: pot, pot_ref
+      character(len=*), intent(in) :: pot, pot_ref
 
       call check_input()
 
@@ -92,7 +93,7 @@ subroutine initialize_integrator(dt, ipimd, inormalmodes, nshake, pot, pot_ref)
       if (.not. associated(mdstep) .and. ipimd /= 3) then
          call fatal_error(__FILE__, __LINE__, 'invalid integrator')
       end if
-   end subroutine
+   end subroutine initialize_integrator
 
    subroutine shiftX(rx, ry, rz, px, py, pz, mass, dt)
       real(DP), intent(inout) :: rx(:, :), ry(:, :), rz(:, :)
@@ -128,7 +129,6 @@ subroutine propagate_nm(x, y, z, px, py, pz, m, dt)
       use mod_const, only: DP, PI
       use mod_general, only: nwalk, natom
       use mod_nhc, only: temp
-      implicit none
       real(DP), intent(inout) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(inout) :: px(:, :), py(:, :), pz(:, :)
       real(DP), intent(in) :: m(:, :)
diff --git a/src/minimizer.F90 b/src/minimizer.F90
index 2ca8104d..15906118 100644
--- a/src/minimizer.F90
+++ b/src/minimizer.F90
@@ -6,6 +6,7 @@
 ! DON'T USE THIS!
 module mod_minimize
    use mod_const, only: DP
+   implicit none
    private
    public :: minimize, gamm, gammthr
    real(DP) :: gamm = 20.D0, gammthr = 1D-10 !minthr=1e-15
@@ -19,7 +20,6 @@ subroutine minimize(x, y, z, fx, fy, fz, eclas)
       use mod_system, only: names, conatom
       use mod_analysis, only: trajout
       use mod_interfaces, only: force_clas
-      implicit none
       real(DP), intent(inout) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :)
       real(DP), intent(inout) :: eclas
diff --git a/src/modules.F90 b/src/modules.F90
index 95e3fd65..4734e23e 100644
--- a/src/modules.F90
+++ b/src/modules.F90
@@ -62,10 +62,10 @@ end subroutine set_natom
 
    subroutine update_simtime(dt)
       use mod_const, only: DP
-      real(DP) :: dt
+      real(DP), intent(in) :: dt
       sim_time = sim_time + dt
    end subroutine update_simtime
-end module
+end module mod_general
 
 ! TODO: Move this to a separate file, and think hard what should be inside this module.
 module mod_system
@@ -91,10 +91,11 @@ subroutine set_atom_names(atnames, num_atom)
 
       allocate (names(num_atom))
       names = atnames
-   end subroutine
+   end subroutine set_atom_names
 end module mod_system
 
 module mod_chars
+   implicit none
    character(len=*), parameter :: CHKNOW = 'If you know what you are doing, &
     &set iknow=1 (namelist general) to proceed.'
 end module mod_chars
diff --git a/src/mpi_wrapper.F90 b/src/mpi_wrapper.F90
index 9a22aa1d..0c58aa53 100644
--- a/src/mpi_wrapper.F90
+++ b/src/mpi_wrapper.F90
@@ -6,6 +6,7 @@ module mod_mpi
 #ifdef USE_MPI
    use mpi
 #endif
+   implicit none
    private
    public :: initialize_mpi, finalize_mpi
    public :: mpi_barrier_wrapper
@@ -78,7 +79,7 @@ subroutine finalize_mpi(error_code)
             call MPI_Finalize(ierr)
          end if
       end if
-   end subroutine
+   end subroutine finalize_mpi
 
    function get_mpi_error_string(mpi_err) result(error_string)
       integer, intent(in) :: mpi_err
diff --git a/src/nosehoover.F90 b/src/nosehoover.F90
index 77b83160..336a461b 100644
--- a/src/nosehoover.F90
+++ b/src/nosehoover.F90
@@ -57,7 +57,7 @@ module mod_nhc
    real(DP) function get_nhcham(natom, nwalk) result(nhcham)
       use mod_system, only: dime
       integer, intent(in) :: natom, nwalk
-      integer iat, inh, iw
+      integer :: iat, inh, iw
 
       nhcham = 0.0D0
       if (imasst == 1) then
@@ -66,13 +66,13 @@ real(DP) function get_nhcham(natom, nwalk) result(nhcham)
             do iw = 1, nwalk
                do iat = 1, natom
                   nhcham = nhcham + &
-                         & pnhx(iat, iw, inh) * pnhx(iat, iw, inh) * 0.5 / Qm(iw) + &
+                         & pnhx(iat, iw, inh) * pnhx(iat, iw, inh) * 0.5D0 / Qm(iw) + &
                          & temp * xi_x(iat, iw, inh)
                   if (dime > 1) nhcham = nhcham + &
-                                       & pnhy(iat, iw, inh) * pnhy(iat, iw, inh) * 0.5 / Qm(iw) + &
+                                       & pnhy(iat, iw, inh) * pnhy(iat, iw, inh) * 0.5D0 / Qm(iw) + &
                                        & temp * xi_y(iat, iw, inh)
                   if (dime > 2) nhcham = nhcham + &
-                                       & pnhz(iat, iw, inh) * pnhz(iat, iw, inh) * 0.5 / Qm(iw) + &
+                                       & pnhz(iat, iw, inh) * pnhz(iat, iw, inh) * 0.5D0 / Qm(iw) + &
                                        & temp * xi_z(iat, iw, inh)
                end do
             end do
@@ -80,14 +80,14 @@ real(DP) function get_nhcham(natom, nwalk) result(nhcham)
 
       else
 
-         iw = 1 !TODO: az bude shake+pimd,tak je tohle treba vyresit
+         iw = 1 !TODO: This needs solving if we ever have SHAKE+PIMD
          do iat = 1, nmolt
             nhcham = nhcham + &
                      & 0.5D0 * pnhx(iat, iw, 1) * pnhx(iat, iw, 1) / ms(iat, 1) + &
                      & dime * natmolt(iat) * temp * xi_x(iat, iw, 1)
             do inh = 2, nchain
                nhcham = nhcham + &
-                        pnhx(iat, iw, inh) * pnhx(iat, iw, inh) * 0.5 / ms(iat, inh) + &
+                        pnhx(iat, iw, inh) * pnhx(iat, iw, inh) * 0.5D0 / ms(iat, inh) + &
                         temp * xi_x(iat, iw, inh)
             end do
          end do
@@ -256,7 +256,7 @@ subroutine set_nhc_massive_masses(nhc_mass)
          end do
 
       end if
-   end subroutine
+   end subroutine set_nhc_massive_masses
 
    subroutine set_nhc_global_masses(mass)
       use mod_system, only: dime
@@ -272,7 +272,7 @@ subroutine set_nhc_global_masses(mass)
             ms(imol, inh) = mass
          end do
       end do
-   end subroutine
+   end subroutine set_nhc_global_masses
 
    subroutine initialize_nhc_momenta(temp)
       use mod_general, only: natom, nwalk
@@ -445,9 +445,10 @@ end function get_nhc_temp
    ! Suzuki-Yoshida split-operator integrator for global NHC
    subroutine shiftNHC_yosh(px, py, pz, amt, dt)
       use mod_system, only: dime
-      real(DP) :: px(:, :), py(:, :), pz(:, :)
-      real(DP) :: amt(:, :), G(MAXCHAIN)
-      real(DP) :: dt, ekin2, AA
+      real(DP), intent(inout) :: px(:, :), py(:, :), pz(:, :)
+      real(DP), intent(in) :: amt(:, :), dt
+      real(DP) :: G(MAXCHAIN)
+      real(DP) :: ekin2, AA
       real(DP) :: wdt, wdt2, wdt4, pscale
       integer :: iw, iat, inh
       integer :: nf, iresp, iyosh
@@ -528,11 +529,10 @@ end subroutine shiftNHC_yosh
    subroutine shiftNHC_yosh_mass(px, py, pz, amt, dt)
       use mod_general
       use mod_shake, only: nshake
-      real(DP) :: px(:, :), py(:, :), pz(:, :)
-      real(DP) :: amt(:, :)
+      real(DP), intent(inout) :: px(:, :), py(:, :), pz(:, :)
+      real(DP), intent(in) :: amt(:, :), dt
       real(DP) :: Gx(MAXCHAIN), Gy(MAXCHAIN), Gz(MAXCHAIN)
-      real(DP) :: dt, AA
-      real(DP) :: wdt, wdt2, wdt4
+      real(DP) :: AA, wdt, wdt2, wdt4
       integer :: iw, iat, inh, istart
       integer :: iresp, iyosh
 
diff --git a/src/pbc.F90 b/src/pbc.F90
index ed027747..dca94a51 100644
--- a/src/pbc.F90
+++ b/src/pbc.F90
@@ -14,7 +14,7 @@ module PBC
 
    subroutine wrap(x, y, z)
       use mod_general, only: nwalk
-      real(DP) :: x(:, :), y(:, :), z(:, :)
+      real(DP), intent(inout) :: x(:, :), y(:, :), z(:, :)
       integer :: i, iat, iat2, iw, iww, iwrap
 
       iwrap = 0
@@ -81,6 +81,6 @@ subroutine wrap(x, y, z)
 
       end do
 
-   end subroutine
+   end subroutine wrap
 
 end module PBC
diff --git a/src/plumed.F90 b/src/plumed.F90
index 229582f8..602a602b 100644
--- a/src/plumed.F90
+++ b/src/plumed.F90
@@ -128,7 +128,6 @@ subroutine force_plumed(x, y, z, fx, fy, fz, eclas)
       use mod_general, only: natom, it, nwalk
       use mod_system, only: am
       use mod_utils, only: c_string
-      implicit none
       real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(in) :: eclas
       real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :)
diff --git a/src/potentials.F90 b/src/potentials.F90
index 9cfa3a26..e4d018b7 100644
--- a/src/potentials.F90
+++ b/src/potentials.F90
@@ -29,7 +29,7 @@ module mod_potentials
       real(DP) :: d0
       real(DP) :: k
       real(DP) :: r0
-   end type
+   end type dw_params
    type(dw_params) :: dw
 
    ! Parameters for two particles bound by Morse potential
@@ -39,21 +39,21 @@ module mod_potentials
       real(DP) :: d0
       real(DP) :: a
       real(DP) :: r0
-   end type
+   end type morse_params
    type(morse_params) :: morse
 
    ! 3D harmonic oscillator
    ! E = 1/2 * (kx * x^2 + ky * y^2 + kz * z^2)
    type :: harm_osc_params
       real(DP) :: kx, ky, kz
-   end type
+   end type harm_osc_params
    type(harm_osc_params) :: ho
 
    ! Two particles bound by harmonic potential
    ! E = 1/2 * k * (r - r0)^2
    type :: harm_rotor_params
       real(DP) :: k, r0
-   end type
+   end type harm_rotor_params
    type(harm_rotor_params) :: hrot
    save
 contains
diff --git a/src/potentials_sh.F90 b/src/potentials_sh.F90
index c9e45c3d..93fa1595 100644
--- a/src/potentials_sh.F90
+++ b/src/potentials_sh.F90
@@ -37,7 +37,7 @@ module mod_potentials_sh
       real(DP) :: a12 = 0.055D0
       real(DP) :: beta12 = 0.6931D0
       real(DP) :: rx = 6.93D0
-   end type
+   end type nai_params
    type(nai_params) :: nai
    save
 contains
diff --git a/src/random.F90 b/src/random.F90
index a16a7a8e..8712ba3e 100644
--- a/src/random.F90
+++ b/src/random.F90
@@ -54,6 +54,7 @@ module mod_random
       use mod_const, only: DP
       use mod_error, only: fatal_error
       use mod_files, only: stdout, stderr
+      implicit none
       private
       public :: gautrg, vranf
       public :: write_prng_state, read_prng_state
@@ -99,7 +100,6 @@ subroutine gautrg(gran, nran, iseed)
 !      subroutines called: vranf, r1mach
 !      m. lewerenz 6/may/90, modified 17/jun/91, mar/95
  
-      implicit none
       integer, intent(in) :: nran
       integer, intent(in), optional :: iseed
       real(DP), intent(inout) :: gran(nran)
@@ -253,7 +253,6 @@ subroutine vranf(ranv, nran, iseed)
 !----------------------------------------------------------------------
 !      subroutines called: xuinit            m. lewerenz may/91 & nov/93
 
-      implicit none
       real(DP), intent(out) :: ranv(nran)
       integer, intent(in) :: nran
       integer, intent(in), optional :: iseed
@@ -351,7 +350,6 @@ subroutine xuinit(y, np, nq, mode, nexec, iseed, init, last)
 !      last  : pointer to the last used number in the table; output
 !      subroutines called : none             m. lewerenz mar/93, mar/98
 
-      implicit none
       real(DP), intent(inout) :: y(np)
       integer, intent(in) :: np, nq, mode, nexec, iseed
       integer, intent(inout) :: init, last
@@ -424,7 +422,6 @@ subroutine xuwarm(y, np, nq, mode, nexec)
 !              random numbers are generated and discarded; input
 !      subroutines called: x                         m. lewerenz mar/98
 
-      implicit none
       ! DHmod, renamed x to y to prevent collision with module
       real(DP), intent(inout) :: y(np)
       real(DP), parameter :: zero = 0.d0, one = 1.d0
@@ -708,7 +705,7 @@ end subroutine random_ints
    ! Taken from:
    ! https://gcc.gnu.org/onlinedocs/gcc-4.9.1/gfortran/RANDOM_005fSEED.html
    integer function lcg(s)
-      integer(INT64) :: s
+      integer(INT64), intent(inout) :: s
 
       if (s == 0) then
          s = 104729
diff --git a/src/remd.F90 b/src/remd.F90
index 0a50aea6..fa9d2db7 100644
--- a/src/remd.F90
+++ b/src/remd.F90
@@ -48,12 +48,15 @@ subroutine remd_swap(x, y, z, px, py, pz, fxc, fyc, fzc, eclas)
       real(DP), intent(inout) :: fxc(:, :), fyc(:, :), fzc(:, :)
       real(DP), intent(inout) :: eclas
       real(DP) :: eclas_new, prob, probs(MAX_REPLICA), ran(MAX_REPLICA), rn
-      integer :: source, ierr, tag_en = 10, tag_swap = 1
+      integer, parameter :: tag_en = 10, tag_swap = 1
+      integer :: source, ierr
       integer :: status(MPI_STATUS_SIZE), i
       integer :: my_rank
-      logical :: lswap = .false., lswaps(MAX_REPLICA)
+      logical :: lswap, lswaps(MAX_REPLICA)
       character(len=100) :: formt
 
+      lswap = .false.
+
       ! Broadcast array of random numbers from rank 0
       my_rank = get_mpi_rank()
       if (my_rank == 0) call vranf(ran, nreplica - 1)
@@ -147,9 +150,9 @@ subroutine swap_replicas(x, y, z, px, py, pz, fxc, fyc, fzc, rank1, rank2)
       integer :: status(MPI_STATUS_SIZE)
       integer :: my_rank
 
-      integer :: tag_x = 11, tag_y = 12, tag_z = 13
-      integer :: tag_px = 114, tag_py = 115, tag_pz = 116
-      integer :: tag_fx = 17, tag_fy = 18, tag_fz = 19
+      integer, parameter :: tag_x = 11, tag_y = 12, tag_z = 13
+      integer, parameter :: tag_px = 114, tag_py = 115, tag_pz = 116
+      integer, parameter :: tag_fx = 17, tag_fy = 18, tag_fz = 19
       integer :: ierr, dest, size1, size2, irank
 
       my_rank = get_mpi_rank()
diff --git a/src/sh_integ.F90 b/src/sh_integ.F90
index 504d1448..c932ae49 100644
--- a/src/sh_integ.F90
+++ b/src/sh_integ.F90
@@ -38,6 +38,7 @@ module mod_sh_integ
    abstract interface
       subroutine sh_integrator(en_array_int, en_array_newint, dotproduct_int, dotproduct_newint, dtp)
          import :: DP
+         implicit none
          ! Interpolated potential energies
          real(DP), intent(in), dimension(:) :: en_array_int, en_array_newint
          ! Interpolated dotproduct = nacm x velocity
diff --git a/src/shake.F90 b/src/shake.F90
index 42df92ce..6e209857 100644
--- a/src/shake.F90
+++ b/src/shake.F90
@@ -25,7 +25,7 @@ subroutine find_hbonds(x, y, z, num_shake)
 
       ! DH fixed threshold for now, we should make it a parameter
       ! Or determine bonds in a better way
-      hbond_len_thr = 2.0 * ANG
+      hbond_len_thr = 2.0_DP * ANG
       ! Here we assume that each H atom has a single bond
       num_shake = count_atoms_by_name(names, 'H', natom)
       ! TODO: Not sure about this, but shake does not
@@ -76,8 +76,8 @@ end subroutine find_hbonds
 
    subroutine shake_init(x, y, z)
       use mod_files, only: stdout
-      real(DP) x(:, :), y(:, :), z(:, :)
-      real(DP) xi, yi, zi, xj, yj, zj
+      real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
+      real(DP) :: xi, yi, zi, xj, yj, zj
       integer :: ixshake, i, j
 
       write (stdout, *) 'Setting distances for SHAKE from XYZ coordinates'
diff --git a/src/surfacehop.F90 b/src/surfacehop.F90
index da4aa731..921c9aeb 100644
--- a/src/surfacehop.F90
+++ b/src/surfacehop.F90
@@ -216,7 +216,7 @@ end subroutine print_sh_input
    subroutine set_current_state(current_state)
       integer, intent(in) :: current_state
       istate = current_state
-   end subroutine
+   end subroutine set_current_state
 
    subroutine check_sh_parameters()
       use mod_utils, only: int_positive, int_nonnegative, int_switch, real_nonnegative
diff --git a/src/tera_mpi_api.F90 b/src/tera_mpi_api.F90
index ae1156a3..ea9fab35 100644
--- a/src/tera_mpi_api.F90
+++ b/src/tera_mpi_api.F90
@@ -47,7 +47,7 @@ module mod_terampi
    subroutine initialize_terachem_interface(tc_server_name)
       use mod_general, only: nwalk, iremd
       use mod_mpi, only: get_mpi_rank
-      character(len=*) :: tc_server_name
+      character(len=*), intent(in) :: tc_server_name
       integer :: i, ierr
 
       if (nwalk > 1) then
diff --git a/src/utils.F90 b/src/utils.F90
index 1336be33..b3e46280 100644
--- a/src/utils.F90
+++ b/src/utils.F90
@@ -28,10 +28,9 @@ end function get_distance
    real(DP) function get_angle(x, y, z, at1, at2, at3, iw) result(angle)
       use mod_const, only: PI
       real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
-      integer, intent(in) :: iw
+      integer, intent(in) :: at1, at2, at3, iw
       real(DP) :: vec1x, vec1y, vec1z
       real(DP) :: vec2x, vec2y, vec2z
-      integer :: at1, at2, at3
 
       if (at1 == at2 .or. at1 == at3 .or. at2 == at3) then
          angle = 0
@@ -53,12 +52,11 @@ end function get_angle
    real(DP) function get_dihedral(x, y, z, at1, at2, at3, at4, iw, shiftdih)
       use mod_const, only: PI
       real(DP), intent(in) :: x(:, :), y(:, :), z(:, :), shiftdih
-      integer, intent(in) :: iw
+      integer, intent(in) :: at1, at2, at3, at4, iw
       real(DP) :: vec1x, vec1y, vec1z
       real(DP) :: vec2x, vec2y, vec2z
       real(DP) :: vec3x, vec3y, vec3z, sign
       real(DP) :: norm1x, norm1y, norm1z, norm2x, norm2y, norm2z
-      integer :: at1, at2, at3, at4
 
       vec1x = x(at1, iw) - x(at2, iw)
       vec1y = y(at1, iw) - y(at2, iw)
@@ -269,7 +267,7 @@ subroutine file_exists_or_exit(fname)
    end subroutine file_exists_or_exit
 
    integer function open_file_for_reading(fname) result(un)
-      character(len=*) :: fname
+      character(len=*), intent(in) :: fname
       character(300) :: errmsg
       integer :: iost
 
@@ -281,7 +279,6 @@ integer function open_file_for_reading(fname) result(un)
    end function open_file_for_reading
 
    real(DP) function ekin_p(px, py, pz, mass, natom, nwalk)
-      implicit none
       real(DP), dimension(:, :), intent(in) :: px, py, pz
       real(DP), dimension(:, :), intent(in) :: mass
       integer, intent(in) :: natom, nwalk
@@ -349,7 +346,7 @@ end subroutine int_switch
    subroutine milisleep(milisec)
       use, intrinsic :: iso_c_binding, only: C_INT, C_INT32_T
       use mod_interfaces, only: usleep
-      integer :: milisec
+      integer, intent(in) :: milisec
       integer(kind=C_INT32_T) :: usec
       integer(kind=C_INT) :: c_err
 
diff --git a/src/vinit.F90 b/src/vinit.F90
index b412f198..79b5f2ac 100644
--- a/src/vinit.F90
+++ b/src/vinit.F90
@@ -37,13 +37,14 @@ module mod_vinit
 ! Adapted (real(DP), no angular velocities, mass)   B. Schmidt, Apr 6, 1995
 !
 !------------------------------------------------------------------------
-   subroutine vinit(TEMP, MASS, vx, vy, vz)
+   subroutine vinit(temp, mass, vx, vy, vz)
       use mod_general, only: natom, nwalk
       use mod_random, only: gautrg
-      real(DP), intent(out) :: vx(:, :), vy(:, :), vz(:, :)
+      real(DP), intent(in) :: temp
       real(DP), intent(in) :: mass(:)
+      real(DP), intent(out) :: vx(:, :), vy(:, :), vz(:, :)
       real(DP) :: rans(3 * size(mass))
-      real(DP) :: TEMP, SIGMA
+      real(DP) :: sigma
       integer :: iw, iat, pom
 
       do iw = 1, nwalk
@@ -264,7 +265,7 @@ subroutine calc_inertia_tensor(x, y, z, masses, iw, I)
       real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(in) :: masses(:)
       integer, intent(in) :: iw ! Bead index
-      real(DP) :: I(0:8)
+      real(DP), intent(out) :: I(0:8)
       real(DP) :: xx, xy, xz, yy, yz, zz
       integer :: iat
 
diff --git a/src/water.F90 b/src/water.F90
index 56ad6157..28a7d3e8 100644
--- a/src/water.F90
+++ b/src/water.F90
@@ -10,7 +10,7 @@ subroutine check_water(natom, names)
       use mod_error, only: fatal_error
       use mod_utils, only: count_atoms_by_name
       integer, intent(in) :: natom
-      character(len=2) :: names(:)
+      character(len=2), intent(in) :: names(:)
       integer :: nH, nO
       integer :: error, iat