diff --git a/sktwocnt/prog/input.f90 b/sktwocnt/prog/input.f90
index 1c636b49..85515dde 100644
--- a/sktwocnt/prog/input.f90
+++ b/sktwocnt/prog/input.f90
@@ -59,8 +59,7 @@ subroutine readInput(inp, fname)
     inp%tCam = .false.
     inp%tGlobalHybrid = .false.
 
-    fp = 14
-    open(fp, file=fname, form="formatted", action="read")
+    open(newunit=fp, file=fname, form="formatted", action="read")
     ! general part
     iLine = 0
 
@@ -348,12 +347,11 @@ subroutine readdata_(fname, cols, data)
     !! auxiliary variable
     integer :: ii
 
-    fp = 12
     iLine = 1
 
     allocate(tmp(maxval(cols)))
 
-    open(fp, file=fname, action="read", form="formatted")
+    open(newunit=fp, file=fname, action="read", form="formatted")
 
     call nextline_(fp, iLine, line)
     read(line, *, iostat=iErr) nGrid
diff --git a/slateratom/lib/output.f90 b/slateratom/lib/output.f90
index c24d7ae5..352906cf 100644
--- a/slateratom/lib/output.f90
+++ b/slateratom/lib/output.f90
@@ -486,8 +486,8 @@ subroutine write_energies_tagged(ekin, enuc, ecoul, exc, econf, etot, zora, eigv
     type(TTaggedwriter) :: twriter
 
     call TTaggedwriter_init(twriter)
-    fp = 95
-    open(fp, file="energies.tag", status="replace", action="write")
+
+    open(newunit=fp, file="energies.tag", status="replace", action="write")
     call writetag(twriter, fp, "zora", zora)
     call writetag(twriter, fp, "kinetic_energy", ekin)
     call writetag(twriter, fp, "nuclear_energy", enuc)
@@ -518,7 +518,7 @@ subroutine write_wave_coeffs_file(max_l, num_alpha, poly_order, cof, alpha, occ,
     real(dp), allocatable :: coeffs(:,:)
 
     call TTaggedwriter_init(twriter)
-    fp = 95
+
     do ll = 0, max_l
       allocate(coeffs(poly_order(ll), num_alpha(ll)))
       do ii = 1, num_alpha(ll) * poly_order(ll)
@@ -526,7 +526,7 @@ subroutine write_wave_coeffs_file(max_l, num_alpha, poly_order, cof, alpha, occ,
           cycle
         end if
         write(fname, "(A,I2.2,A,A)") "coeffs_", ii + ll, orbnames(ll), ".tag"
-        open(fp, file=fname, status="replace", action="write")
+        open(newunit=fp, file=fname, status="replace", action="write")
         call writetag(twriter, fp, "exponents", alpha(ll, :num_alpha(ll)))
         call convcoeffs(cof(1, ll, :, ii), alpha(ll, :num_alpha(ll)), ll, coeffs)
         call writetag(twriter, fp, "coefficients", coeffs)