diff --git a/src/CLI.f90 b/src/CLI.f90
index cc73365d6..ee1555cb2 100644
--- a/src/CLI.f90
+++ b/src/CLI.f90
@@ -18,7 +18,7 @@ module CLI
   use parallelization
   use system_routines
 
-  implicit none
+  implicit none(type,external)
   private
   integer,                       public, protected :: &
     CLI_restartInc = 0                                                                              !< Increment at which calculation starts
diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90
index c905164c2..ee1e3cf8c 100644
--- a/src/HDF5_utilities.f90
+++ b/src/HDF5_utilities.f90
@@ -18,7 +18,11 @@ module HDF5_utilities
   use prec
   use parallelization
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
 !--------------------------------------------------------------------------------------------------
diff --git a/src/IO.f90 b/src/IO.f90
index 9b9b9d563..31bb0ed2d 100644
--- a/src/IO.f90
+++ b/src/IO.f90
@@ -12,7 +12,7 @@ module IO
 
   use prec
 
-  implicit none
+  implicit none(type,external)
   private
 
   character(len=*), parameter, public :: &
diff --git a/src/LAPACK_interface.f90 b/src/LAPACK_interface.f90
index 43ddb5771..cc451b59c 100644
--- a/src/LAPACK_interface.f90
+++ b/src/LAPACK_interface.f90
@@ -8,7 +8,7 @@ module LAPACK_interface
 
     pure subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)
       use prec
-      implicit none
+      implicit none(type,external)
 
       character,   intent(in)                             :: jobvl,jobvr
       integer,     intent(in)                             :: n,lda,ldvl,ldvr,lwork
@@ -22,7 +22,7 @@ module LAPACK_interface
 
     pure subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
       use prec
-      implicit none
+      implicit none(type,external)
 
       integer,     intent(in)                             :: n,nrhs,lda,ldb
       real(pReal), intent(inout), dimension(lda,n)        :: a
@@ -33,7 +33,7 @@ module LAPACK_interface
 
     pure subroutine dgetrf(m,n,a,lda,ipiv,info)
       use prec
-      implicit none
+      implicit none(type,external)
 
       integer,     intent(in)                             :: m,n,lda
       real(pReal), intent(inout), dimension(lda,n)        :: a
@@ -43,7 +43,7 @@ module LAPACK_interface
 
     pure subroutine dgetri(n,a,lda,ipiv,work,lwork,info)
       use prec
-      implicit none
+      implicit none(type,external)
 
       integer,     intent(in)                             :: n,lda,lwork
       real(pReal), intent(inout), dimension(lda,n)        :: a
@@ -54,7 +54,7 @@ module LAPACK_interface
 
     pure subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info)
       use prec
-      implicit none
+      implicit none(type,external)
 
       character,   intent(in)                             :: jobz,uplo
       integer,     intent(in)                             :: n,lda,lwork
diff --git a/src/Marc/DAMASK_Marc.f90 b/src/Marc/DAMASK_Marc.f90
index 3cfaf7741..c1f5b1d68 100644
--- a/src/Marc/DAMASK_Marc.f90
+++ b/src/Marc/DAMASK_Marc.f90
@@ -25,7 +25,7 @@ module DAMASK_interface
   use ifport, only: &
     CHDIR
 
-  implicit none
+  implicit none(type,external)
   private
 
   logical,          protected, public :: symmetricSolver
@@ -210,7 +210,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
   use materialpoint_Marc
   use OMP_LIB
 
-  implicit none
+  implicit none(type,external)
   integer(pI64),                         intent(in) :: &                                            ! according to MSC.Marc 2012 Manual D
     ngens, &                                                                                        !< size of stress-strain law
     nn, &                                                                                           !< integration point number
@@ -382,7 +382,7 @@ subroutine flux(f,ts,n,time)
   use homogenization
   use discretization_Marc
 
-  implicit none
+  implicit none(type,external)
   real(pReal),   dimension(6),  intent(in) :: &
     ts
   integer(pI64), dimension(10), intent(in) :: &
@@ -410,7 +410,7 @@ subroutine uedinc(inc,incsub)
   use materialpoint_Marc
   use discretization_Marc
 
-  implicit none
+  implicit none(type,external)
   integer(pI64), intent(in) :: inc, incsub
 
   integer :: n, nqncomp, nqdatatype
diff --git a/src/Marc/discretization_Marc.f90 b/src/Marc/discretization_Marc.f90
index 8886e3ebb..a0667d049 100644
--- a/src/Marc/discretization_Marc.f90
+++ b/src/Marc/discretization_Marc.f90
@@ -17,7 +17,7 @@ module discretization_Marc
   use geometry_plastic_nonlocal
   use results
 
-  implicit none
+  implicit none(type,external)
   private
 
   real(pReal),                         public, protected :: &
diff --git a/src/Marc/element.f90 b/src/Marc/element.f90
index b53322f34..2940704bd 100644
--- a/src/Marc/element.f90
+++ b/src/Marc/element.f90
@@ -5,7 +5,7 @@
 module element
   use IO
 
-  implicit none
+  implicit none(type,external)
   private
 
 !---------------------------------------------------------------------------------------------------
diff --git a/src/Marc/materialpoint_Marc.f90 b/src/Marc/materialpoint_Marc.f90
index 2ad98109d..72c718a87 100644
--- a/src/Marc/materialpoint_Marc.f90
+++ b/src/Marc/materialpoint_Marc.f90
@@ -23,7 +23,7 @@ module materialpoint_Marc
   use discretization
   use discretization_Marc
 
-  implicit none
+  implicit none(type,external)
   private
 
   real(pReal), dimension (:,:,:),   allocatable, private :: &
diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90
index 67abfdd90..e7d64770d 100644
--- a/src/YAML_parse.f90
+++ b/src/YAML_parse.f90
@@ -12,7 +12,7 @@ module YAML_parse
   use system_routines
 #endif
 
-  implicit none
+  implicit none(type,external)
   private
 
   public :: &
@@ -24,7 +24,7 @@ module YAML_parse
 
     subroutine to_flow_C(flow,length_flow,mixed) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR, C_PTR
-      implicit none
+      implicit none(type,external)
 
       type(C_PTR), intent(out) :: flow
       integer(C_INT), intent(out) :: length_flow
diff --git a/src/YAML_types.f90 b/src/YAML_types.f90
index e685209e6..e9bdff88d 100644
--- a/src/YAML_types.f90
+++ b/src/YAML_types.f90
@@ -11,7 +11,7 @@ module YAML_types
   use IO
   use prec
 
-  implicit none
+  implicit none(type,external)
   private
 
   type, abstract, public :: tNode
diff --git a/src/config.f90 b/src/config.f90
index 78723c19f..5c31e3b25 100644
--- a/src/config.f90
+++ b/src/config.f90
@@ -9,7 +9,7 @@ module config
   use results
   use parallelization
 
-  implicit none
+  implicit none(type,external)
   private
 
   class(tNode), pointer, public :: &
diff --git a/src/constants.f90 b/src/constants.f90
index 1effb8d95..7d1892e9d 100644
--- a/src/constants.f90
+++ b/src/constants.f90
@@ -5,7 +5,7 @@
 module constants
   use prec
 
-  implicit none
+  implicit none(type,external)
   public
 
   real(pReal), parameter :: &
diff --git a/src/discretization.f90 b/src/discretization.f90
index 95fd39fd5..f24b0eadf 100644
--- a/src/discretization.f90
+++ b/src/discretization.f90
@@ -7,7 +7,7 @@ module discretization
   use prec
   use results
 
-  implicit none
+  implicit none(type,external)
   private
 
   integer,     public, protected :: &
diff --git a/src/geometry_plastic_nonlocal.f90 b/src/geometry_plastic_nonlocal.f90
index 571ca2e6a..09c40f8b3 100644
--- a/src/geometry_plastic_nonlocal.f90
+++ b/src/geometry_plastic_nonlocal.f90
@@ -9,7 +9,7 @@ module geometry_plastic_nonlocal
   use prec
   use results
 
-  implicit none
+  implicit none(type,external)
   public
 
   integer, protected :: &
diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90
index 0f16d33aa..4395c1581 100644
--- a/src/grid/DAMASK_grid.f90
+++ b/src/grid/DAMASK_grid.f90
@@ -30,7 +30,11 @@ program DAMASK_grid
   use grid_thermal_spectral
   use results
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
 
   type :: tLoadCase
     type(tRotation)          :: rot                                                                !< rotation of BC
diff --git a/src/grid/FFTW.f90 b/src/grid/FFTW.f90
index e552e78e0..5df432061 100644
--- a/src/grid/FFTW.f90
+++ b/src/grid/FFTW.f90
@@ -5,7 +5,7 @@
 module FFTW3
   use, intrinsic :: ISO_C_binding
 
-  implicit none
+  implicit none(type,external)
   public
 
   include 'fftw3-mpi.f03'
diff --git a/src/grid/VTI.f90 b/src/grid/VTI.f90
index 53be9f022..72f57fb86 100644
--- a/src/grid/VTI.f90
+++ b/src/grid/VTI.f90
@@ -8,7 +8,7 @@ module VTI
   use base64
   use IO
 
-  implicit none
+  implicit none(type,external)
   private
 
   public :: &
diff --git a/src/grid/base64.f90 b/src/grid/base64.f90
index 448f57d2a..6e580f043 100644
--- a/src/grid/base64.f90
+++ b/src/grid/base64.f90
@@ -7,7 +7,7 @@ module base64
   use prec
   use IO
 
-  implicit none
+  implicit none(type,external)
   private
 
   character(len=*), parameter :: &
diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90
index 2bac30f0a..158ee0a8d 100644
--- a/src/grid/discretization_grid.f90
+++ b/src/grid/discretization_grid.f90
@@ -23,7 +23,11 @@ module discretization_grid
   use discretization
   use geometry_plastic_nonlocal
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   integer,     dimension(3), public, protected :: &
diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90
index dc21a1a2c..8b544ec87 100644
--- a/src/grid/grid_damage_spectral.f90
+++ b/src/grid/grid_damage_spectral.f90
@@ -22,7 +22,11 @@ module grid_damage_spectral
   use YAML_types
   use config
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   type :: tNumerics
diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90
index 405b10845..4bb705807 100644
--- a/src/grid/grid_mech_FEM.f90
+++ b/src/grid/grid_mech_FEM.f90
@@ -27,7 +27,12 @@ module grid_mechanical_FEM
   use discretization
   use discretization_grid
 
+
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   type(tSolutionParams) :: params
diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90
index 7931be0c7..dec677df8 100644
--- a/src/grid/grid_mech_spectral_basic.f90
+++ b/src/grid/grid_mech_spectral_basic.f90
@@ -26,7 +26,11 @@ module grid_mechanical_spectral_basic
   use homogenization
   use discretization_grid
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   type(tSolutionParams) :: params
diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90
index 198e333d5..d0d229c46 100644
--- a/src/grid/grid_mech_spectral_polarisation.f90
+++ b/src/grid/grid_mech_spectral_polarisation.f90
@@ -26,7 +26,11 @@ module grid_mechanical_spectral_polarisation
   use homogenization
   use discretization_grid
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   type(tSolutionParams) :: params
diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90
index 8483fe336..ff0c5b403 100644
--- a/src/grid/grid_thermal_spectral.f90
+++ b/src/grid/grid_thermal_spectral.f90
@@ -25,7 +25,11 @@ module grid_thermal_spectral
   use YAML_types
   use config
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   type :: tNumerics
diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90
index 1c8842a60..95e75b854 100644
--- a/src/grid/spectral_utilities.f90
+++ b/src/grid/spectral_utilities.f90
@@ -22,7 +22,12 @@ module spectral_utilities
   use discretization
   use homogenization
 
+
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
 !--------------------------------------------------------------------------------------------------
diff --git a/src/grid/zlib.f90 b/src/grid/zlib.f90
index 10a62a55b..f0fc274ed 100644
--- a/src/grid/zlib.f90
+++ b/src/grid/zlib.f90
@@ -5,7 +5,7 @@
 module zlib
   use prec
 
-  implicit none
+  implicit none(type,external)
   private
 
   public :: &
@@ -15,7 +15,7 @@ module zlib
 
     subroutine inflate_C(s_deflated,s_inflated,deflated,inflated) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_SIGNED_CHAR, C_INT64_T
-      implicit none
+      implicit none(type,external)
 
       integer(C_INT64_T),                            intent(in)  :: s_deflated,s_inflated
       integer(C_SIGNED_CHAR), dimension(s_deflated), intent(in)  :: deflated
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 290f3de9b..b1e462794 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -18,7 +18,7 @@ module homogenization
   use results
   use lattice
 
-  implicit none
+  implicit none(type,external)
   private
 
   type :: tState
diff --git a/src/lattice.f90 b/src/lattice.f90
index 800b33dc3..7bd846d0f 100644
--- a/src/lattice.f90
+++ b/src/lattice.f90
@@ -13,7 +13,7 @@ module lattice
   use math
   use rotations
 
-  implicit none
+  implicit none(type,external)
   private
 
 !--------------------------------------------------------------------------------------------------
diff --git a/src/material.f90 b/src/material.f90
index eafc411cd..1ec844382 100644
--- a/src/material.f90
+++ b/src/material.f90
@@ -14,7 +14,7 @@ module material
   use discretization
   use YAML_types
 
-  implicit none
+  implicit none(type,external)
   private
 
   type, public :: tRotationContainer
diff --git a/src/materialpoint.f90 b/src/materialpoint.f90
index d83307738..5b8b690e0 100644
--- a/src/materialpoint.f90
+++ b/src/materialpoint.f90
@@ -31,7 +31,7 @@ module materialpoint
   use discretization_grid
 #endif
 
-  implicit none
+  implicit none(type,external)
   public
 
 contains
diff --git a/src/math.f90 b/src/math.f90
index 45d0da972..07b657cab 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -12,7 +12,7 @@ module math
   use YAML_types
   use LAPACK_interface
 
-  implicit none
+  implicit none(type,external)
   public
 #if __INTEL_COMPILER >= 1900
   ! do not make use of associated entities available to other modules
diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90
index f0f2a872e..5489ac36e 100644
--- a/src/mesh/DAMASK_mesh.f90
+++ b/src/mesh/DAMASK_mesh.f90
@@ -20,7 +20,7 @@ program DAMASK_mesh
   use FEM_Utilities
   use mesh_mechanical_FEM
 
-  implicit none
+  implicit none(type,external)
 
   type :: tLoadCase
     real(pReal)  :: time                   = 0.0_pReal                                              !< length of increment
diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90
index f519dc377..a2217847a 100644
--- a/src/mesh/FEM_quadrature.f90
+++ b/src/mesh/FEM_quadrature.f90
@@ -5,7 +5,7 @@
 module FEM_quadrature
   use prec
 
-  implicit none
+  implicit none(type,external)
   private
 
   integer, parameter :: &
diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90
index 0b633bef7..c97c2bbe2 100644
--- a/src/mesh/FEM_utilities.f90
+++ b/src/mesh/FEM_utilities.f90
@@ -21,7 +21,11 @@ module FEM_utilities
   use homogenization
   use FEM_quadrature
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   logical,     public             :: cutBack = .false.                                              !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
@@ -65,6 +69,11 @@ module FEM_utilities
     type(tComponentBC), allocatable, dimension(:) :: componentBC
   end type tFieldBC
 
+  external :: &                                                                                     ! ToDo: write interfaces
+    PetscSectionGetFieldComponents, &
+    PetscSectionGetFieldDof, &
+    PetscSectionGetFieldOffset
+
   public :: &
     FEM_utilities_init, &
     utilities_constitutiveResponse, &
diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90
index 206d6b5bc..a53c50655 100644
--- a/src/mesh/discretization_mesh.f90
+++ b/src/mesh/discretization_mesh.f90
@@ -25,7 +25,11 @@ module discretization_mesh
   use YAML_types
   use prec
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   integer, public, protected :: &
@@ -52,6 +56,11 @@ module discretization_mesh
   real(pReal), dimension(:,:,:), allocatable :: &
     mesh_ipCoordinates                                                                              !< IP x,y,z coordinates (after deformation!)
 
+  external :: &
+#ifdef PETSC_USE_64BIT_INDICES
+    DMDestroy, &
+#endif
+    DMView                                                                                          ! ToDo: write interface
   public :: &
     discretization_mesh_init, &
     mesh_FEM_build_ipVolumes, &
@@ -242,10 +251,10 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
     call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc)
     CHKERRQ(err_PETSc)
     qOffset = 0
-    do qPt = 1, mesh_maxNips
-      do dirI = 1, dimPlex
+    do qPt = 1_pPETSCINT, mesh_maxNips
+      do dirI = 1_pPETSCINT, dimPlex
         mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI)
-        do dirJ = 1, dimPlex
+        do dirJ = 1_pPETSCINT, dimPlex
           mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + &
                                                 pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0_pReal)
         enddo
diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90
index 4ba3f2690..3b162d97b 100644
--- a/src/mesh/mesh_mech_FEM.f90
+++ b/src/mesh/mesh_mech_FEM.f90
@@ -26,7 +26,11 @@ module mesh_mechanical_FEM
   use homogenization
   use math
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
 !--------------------------------------------------------------------------------------------------
@@ -67,6 +71,18 @@ module mesh_mechanical_FEM
   logical :: ForwardData
   real(pReal), parameter :: eps = 1.0e-18_pReal
 
+  external :: &                                                                                     ! ToDo: write interfaces
+#ifdef PETSC_USE_64BIT_INDICES
+    ISDestroy, &
+#endif
+    PetscSectionGetNumFields, &
+    PetscFESetQuadrature, &
+    PetscFEGetDimension, &
+    PetscFEDestroy, &
+    PetscSectionGetDof, &
+    PetscFEGetDualSpace, &
+    PetscDualSpaceGetFunctional
+
   public :: &
     FEM_mechanical_init, &
     FEM_mechanical_solution, &
@@ -359,7 +375,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
   CHKERRQ(err_PETSc)
   call VecWAXPY(x_local,1.0_pReal,xx_local,solution_local,err_PETSc)
   CHKERRQ(err_PETSc)
-  do field = 1, dimPlex; do face = 1, mesh_Nboundaries
+  do field = 1_pPETSCINT, dimPlex; do face = 1_pPETSCINT, mesh_Nboundaries
     if (params%fieldBC%componentBC(field)%Mask(face)) then
       call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
       if (bcSize > 0) then
@@ -374,7 +390,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
 
 !--------------------------------------------------------------------------------------------------
 ! evaluate field derivatives
-  do cell = cellStart, cellEnd-1                                                                    !< loop over all elements
+  do cell = cellStart, cellEnd-1_pPETSCINT                                                          !< loop over all elements
 
     call PetscSectionGetNumFields(section,numFields,err_PETSc)
     CHKERRQ(err_PETSc)
@@ -383,11 +399,11 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
     call  DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc)
     CHKERRQ(err_PETSc)
     IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
-    do qPt = 0, nQuadrature-1
-      m = cell*nQuadrature + qPt+1
+    do qPt = 0_pPETSCINT, nQuadrature-1_pPETSCINT
+      m = cell*nQuadrature + qPt+1_pPETSCINT
       BMat = 0.0_pReal
-      do basis = 0, nBasis-1
-        do comp = 0, dimPlex-1
+      do basis = 0_pPETSCINT, nBasis-1_pPETSCINT
+        do comp = 0_pPETSCINT, dimPlex-1_pPETSCINT
           cidx = basis*dimPlex+comp
           i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp
           BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = &
@@ -425,11 +441,11 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
     CHKERRQ(err_PETSc)
     IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
     f_scal = 0.0_pReal
-    do qPt = 0, nQuadrature-1
-      m = cell*nQuadrature + qPt+1
+    do qPt = 0_pPETSCINT, nQuadrature-1_pPETSCINT
+      m = cell*nQuadrature + qPt+1_pPETSCINT
       BMat = 0.0_pReal
-      do basis = 0, nBasis-1
-        do comp = 0, dimPlex-1
+      do basis = 0_pPETSCINT, nBasis-1_pPETSCINT
+        do comp = 0_pPETSCINT, dimPlex-1_pPETSCINT
           cidx = basis*dimPlex+comp
           i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp
           BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = &
@@ -531,11 +547,11 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
     MatB = 0.0_pReal
     FAvg = 0.0_pReal
     BMatAvg = 0.0_pReal
-    do qPt = 0, nQuadrature-1
-      m = cell*nQuadrature + qPt + 1
+    do qPt = 0_pPETSCINT, nQuadrature-1_pPETSCINT
+      m = cell*nQuadrature + qPt + 1_pPETSCINT
       BMat = 0.0_pReal
-      do basis = 0, nBasis-1
-        do comp = 0, dimPlex-1
+      do basis = 0_pPETSCINT, nBasis-1_pPETSCINT
+        do comp = 0_pPETSCINT, dimPlex-1_pPETSCINT
           cidx = basis*dimPlex+comp
           i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp
           BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = &
@@ -741,7 +757,7 @@ subroutine FEM_mechanical_updateCoords()
   call PetscDSGetTabulation(mechQuad,0_pPETSCINT,basisField,basisFieldDer,err_PETSc)
   CHKERRQ(err_PETSc)
   allocate(ipCoords(3,nQuadrature,mesh_NcpElems),source=0.0_pReal)
-  do c=cellStart,cellEnd-1
+  do c=cellStart,cellEnd-1_pPETSCINT
     qOffset=0
     call DMPlexVecGetClosure(dm_local,section,x_local,c,x_scal,err_PETSc)                           !< get nodal coordinates of each element
     CHKERRQ(err_PETSc)
diff --git a/src/parallelization.f90 b/src/parallelization.f90
index d631711d8..5c7e9bff3 100644
--- a/src/parallelization.f90
+++ b/src/parallelization.f90
@@ -18,7 +18,11 @@ module parallelization
 
   use prec
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
 #ifndef PETSC
diff --git a/src/phase.f90 b/src/phase.f90
index 1403fef22..f304707b9 100644
--- a/src/phase.f90
+++ b/src/phase.f90
@@ -19,7 +19,7 @@ module phase
   use HDF5
   use HDF5_utilities
 
-  implicit none
+  implicit none(type,external)
   private
 
   type :: tState
diff --git a/src/polynomials.f90 b/src/polynomials.f90
index e04db16f5..eeb4152f2 100644
--- a/src/polynomials.f90
+++ b/src/polynomials.f90
@@ -8,7 +8,7 @@ module polynomials
   use YAML_parse
   use YAML_types
 
-  implicit none
+  implicit none(type,external)
   private
 
   type, public :: tPolynomial
diff --git a/src/prec.f90 b/src/prec.f90
index 73563d4f2..963391c39 100644
--- a/src/prec.f90
+++ b/src/prec.f90
@@ -15,7 +15,7 @@ module prec
   use PETScSys
 #endif
 
-  implicit none
+  implicit none(type,external)
   public
 
   ! https://stevelionel.com/drfortran/2017/03/27/doctor-fortran-in-it-takes-all-kinds
diff --git a/src/quit.f90 b/src/quit.f90
index aa8bcb0fb..4fc6bf11b 100644
--- a/src/quit.f90
+++ b/src/quit.f90
@@ -12,8 +12,14 @@ subroutine quit(stop_id)
 #endif
   use HDF5
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
+
   integer, intent(in) :: stop_id
+
   integer, dimension(8) :: dateAndTime
   integer :: err_HDF5
   integer(MPI_INTEGER_KIND) :: err_MPI
diff --git a/src/results.f90 b/src/results.f90
index cf1b4de5c..17efc87e3 100644
--- a/src/results.f90
+++ b/src/results.f90
@@ -21,7 +21,11 @@ module results
   use DAMASK_interface
 #endif
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   integer(HID_T) :: resultsFile
diff --git a/src/rotations.f90 b/src/rotations.f90
index 008078c15..5c8677c81 100644
--- a/src/rotations.f90
+++ b/src/rotations.f90
@@ -50,7 +50,7 @@ module rotations
   use IO
   use math
 
-  implicit none
+  implicit none(type,external)
   private
 
   real(pReal), parameter :: P = -1.0_pReal                                                          !< parameter for orientation conversion.
@@ -796,7 +796,7 @@ end function conjugateQuaternion
 subroutine selfTest()
 
   type(tRotation)                 :: R
-  real(pReal), dimension(4)       :: qu, ax
+  real(pReal), dimension(4)       :: qu
   real(pReal), dimension(3)       :: x, eu, v3
   real(pReal), dimension(3,3)     :: om, t33
   real(pReal), dimension(3,3,3,3) :: t3333
diff --git a/src/signals.f90 b/src/signals.f90
index 7db101f94..3f0397d3d 100644
--- a/src/signals.f90
+++ b/src/signals.f90
@@ -6,7 +6,7 @@ module signals
   use prec
   use system_routines
 
-  implicit none
+  implicit none(type,external)
   private
 
   logical, volatile, public, protected :: &
diff --git a/src/system_routines.f90 b/src/system_routines.f90
index a74c2e417..3ce6ba6ce 100644
--- a/src/system_routines.f90
+++ b/src/system_routines.f90
@@ -7,7 +7,7 @@ module system_routines
 
   use prec
 
-  implicit none
+  implicit none(type,external)
   private
 
   public :: &
@@ -27,7 +27,7 @@ module system_routines
 
     function setCWD_C(cwd) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
-      implicit none
+      implicit none(type,external)
 
       integer(C_INT) :: setCWD_C
       character(kind=C_CHAR), dimension(*), intent(in) :: cwd
@@ -36,7 +36,7 @@ module system_routines
     subroutine getCWD_C(cwd, stat) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
       use prec
-      implicit none
+      implicit none(type,external)
 
       character(kind=C_CHAR), dimension(pPathLen+1), intent(out) :: cwd                             ! NULL-terminated array
       integer(C_INT),                                intent(out) :: stat
@@ -45,7 +45,7 @@ module system_routines
     subroutine getHostName_C(hostname, stat) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
       use prec
-      implicit none
+      implicit none(type,external)
 
       character(kind=C_CHAR), dimension(pStringLen+1), intent(out) :: hostname                        ! NULL-terminated array
       integer(C_INT),                                  intent(out) :: stat
@@ -54,7 +54,7 @@ module system_routines
     subroutine getUserName_C(username, stat) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
       use prec
-      implicit none
+      implicit none(type,external)
 
       character(kind=C_CHAR), dimension(pStringLen+1), intent(out) :: username                        ! NULL-terminated array
       integer(C_INT),                                  intent(out) :: stat
@@ -62,28 +62,28 @@ module system_routines
 
     subroutine signalint_C(handler) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_FUNPTR
-      implicit none
+      implicit none(type,external)
 
       type(C_FUNPTR), intent(in), value :: handler
     end subroutine signalint_C
 
     subroutine signalusr1_C(handler) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_FUNPTR
-      implicit none
+      implicit none(type,external)
 
       type(C_FUNPTR), intent(in), value :: handler
     end subroutine signalusr1_C
 
     subroutine signalusr2_C(handler) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_FUNPTR
-      implicit none
+      implicit none(type,external)
 
       type(C_FUNPTR), intent(in), value :: handler
     end subroutine signalusr2_C
 
     subroutine free_C(ptr) bind(C,name='free')
       use, intrinsic :: ISO_C_Binding, only: C_PTR
-      implicit none
+      implicit none(type,external)
 
       type(C_PTR), value :: ptr
     end subroutine free_C