diff --git a/LICENSE b/LICENSE
index 15e905f..b2713e9 100644
--- a/LICENSE
+++ b/LICENSE
@@ -19,3 +19,13 @@ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 SOFTWARE.
+
+
+## License for src/modspainv.f90
+
+The code in modspainv.f90 is based on code originally written by Karin Meyer
+(didgeridoo.une.edu.au/womwiki/doku.php?id=fortran:fortran). The incorporated
+and modified code is re-licensed under the MIT license with express permission
+from Karin Meyer. We are grateful for her permission to use and adapt this
+work.
+
diff --git a/README.md b/README.md
index b7e198b..fda14b8 100644
--- a/README.md
+++ b/README.md
@@ -37,6 +37,12 @@ The brief documentation is available in the directory *doc* (see *mainpage.md*).
 
 
 ## Acknowledgements  
+
+The code in modspainv.f90 is based on code originally written by Karin Meyer
+(didgeridoo.une.edu.au/womwiki/doku.php?id=fortran:fortran). The incorporated and
+modified code is re-licensed under the MIT license with express permission from
+Karin Meyer. We are grateful for her permission to use and adapt this work.
+
 This library was inspired by several sources:  
 
 
diff --git a/src/modspainv.f90 b/src/modspainv.f90
index f010fe0..19387c0 100644
--- a/src/modspainv.f90
+++ b/src/modspainv.f90
@@ -6,7 +6,11 @@
 !Support of sparse inverse of SPSD matrices by implementing the S. D. Kachman modifications
 !(https://www.ars.usda.gov/ARSUserFiles/80420530/MTDFREML/MTDFMan.pdf ; Chapter 6)
 
-!Rewritten for my purposes
+!The code in modspainv.f90 is based on code originally written by Karin Meyer
+!(didgeridoo.une.edu.au/womwiki/doku.php?id=fortran:fortran). The incorporated and 
+!modified code is re-licensed under the MIT license with express permission from
+!Karin Meyer. We are grateful for her permission to use and adapt this work.
+
 
 module modspainv
 #if (_DP==0)
diff --git a/src/modsparse_inv_int.f90 b/src/modsparse_inv_int.f90
index 08ad10b..9de1e0f 100644
--- a/src/modsparse_inv_int.f90
+++ b/src/modsparse_inv_int.f90
@@ -92,40 +92,40 @@ subroutine gsfct(neqns, xlnz, lnz, xnzsub, nzsub, diag, iflag)
   allocate (link(neqns), source=0)
   allocate (temp(neqns), source=0.0_wp)
   iflag = 0
-!c        --------------------------------------------
-!c        compute column l(*,j) for j = 1,...., neqns.
-!c        --------------------------------------------
+  !--------------------------------------------
+  !compute column l(*,j) for j = 1,...., neqns.
+  !--------------------------------------------
   do j = 1, neqns
-!c           -------------------------------------------
-!c           for each column l(*,k) that affects l(*,j).
-!c           -------------------------------------------
+   !-------------------------------------------
+   !for each column l(*,k) that affects l(*,j).
+   !-------------------------------------------
    diagj = 0.0d0
    newk = link(j)
    k = newk
    do while (k /= 0)
     newk = link(k)
-!c              ---------------------------------------
-!c              outer product modification of l(*,j) by
-!c              l(*,k) starting at first(k) of l(*,k).
-!c              ---------------------------------------
+    !---------------------------------------
+    !outer product modification of l(*,j) by
+    !l(*,k) starting at first(k) of l(*,k).
+    !---------------------------------------
     kfirst = first(k)
     ljk = lnz(kfirst)
     diagj = diagj + ljk*ljk
     istrt = kfirst + 1
     istop = xlnz(k + 1) - 1
     if (istop >= istrt) then
-!c                 ------------------------------------------
-!c                 before modification, update vectors first,
-!c                 and link for future modification steps.
-!c                 ------------------------------------------
+     !------------------------------------------
+     !before modification, update vectors first,
+     !and link for future modification steps.
+     !------------------------------------------
      first(k) = istrt
      i = xnzsub(k) + (kfirst - xlnz(k)) + 1
      isub = nzsub(i)
      link(k) = link(isub)
      link(isub) = k
-!c                 ---------------------------------------
-!c                 the actual mod is saved in vector temp.
-!c                 ---------------------------------------
+     !---------------------------------------
+     !the actual mod is saved in vector temp.
+     !---------------------------------------
      do ii = istrt, istop
       isub = nzsub(i)
       temp(isub) = temp(isub) + lnz(ii)*ljk
@@ -134,23 +134,23 @@ subroutine gsfct(neqns, xlnz, lnz, xnzsub, nzsub, diag, iflag)
     end if
     k = newk
    end do
-!c           ----------------------------------------------
-!c           apply the modifications accumulated in temp to
-!c           column l(*,j).
-!c           ----------------------------------------------
+   !----------------------------------------------
+   !apply the modifications accumulated in temp to
+   !column l(*,j).
+   !----------------------------------------------
    diagj = diag(j) - diagj
    if (diag(j) < tol) then
-!c              ------------------------------------------------------
-!c              error - zero diagonal element
-!c              ------------------------------------------------------
+    !------------------------------------------------------
+    !error - zero diagonal element
+    !------------------------------------------------------
     diagj = 0.d0
     diag(j) = 0.d0
     iflag = iflag + 1
    else
     if (diagj < tol*diag(j)) then
-!c              ------------------------------------------------------
-!c              error - zero or negative square root in factorization.
-!c              ------------------------------------------------------
+     !------------------------------------------------------
+     !error - zero or negative square root in factorization.
+     !------------------------------------------------------
      diagj = 0.d0
      diag(j) = 0.d0
      iflag = iflag + 1