diff --git a/src/m_common/upwinv.f90 b/src/m_common/upwinv.f90
index b3abd1f..457529f 100644
--- a/src/m_common/upwinv.f90
+++ b/src/m_common/upwinv.f90
@@ -16,8 +16,15 @@
 subroutine upwinv(nel, jel, indt, nelorb, winv, v, psi)
     use constants, only: yes_ontarget
     implicit none
-    integer nel, indt, i, j, jel, nelorb
-    real*8 winv(nel, indt), psi(indt, nel), v(nel)
+
+    ! argument parameters
+    integer, intent(in) :: nel, indt, jel, nelorb
+    real*8, intent(in) :: psi(indt, nel), v(nel)
+    real*8, intent(inout) :: winv(nel, indt)
+
+    ! local variables
+    integer i, j
+
 #ifdef _OFFLOAD
 !$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
 #endif
@@ -43,11 +50,18 @@ subroutine upwinv(nel, jel, indt, nelorb, winv, v, psi)
     end do
     return
 end
+
 subroutine upwinv_complex(nel, jel, indt, nelorb, winv, v, psi)
     use constants, only: yes_ontarget
     implicit none
-    integer nel, indt, i, j, jel, nelorb
-    complex*16 winv(nel, indt), psi(indt, nel), v(nel)
+
+    ! argument parameters
+    integer, intent(in) :: nel, indt, jel, nelorb
+    complex*16, intent(in) :: psi(indt, nel), v(nel)
+    complex*16, intent(inout) :: winv(nel, indt)
+
+    ! local variables
+    integer i, j
 #ifdef _OFFLOAD
 !$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
 #endif
diff --git a/src/m_common/upwinvp.f90 b/src/m_common/upwinvp.f90
index 3062009..fe835f2 100644
--- a/src/m_common/upwinvp.f90
+++ b/src/m_common/upwinvp.f90
@@ -13,11 +13,19 @@
 ! You should have received a copy of the GNU General Public License
 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
 
+!> subroutine to update the matrix winv, for real(8)
 subroutine upwinvp(nel, indt, winv, ainv, ainvn, psi)
     use constants, only: yes_ontarget
     implicit none
-    integer nel, indt, i, j
-    real*8 winv(nel, indt), psi(indt, nel, 2), ainv(nel), ainvn(nel)
+
+    ! argument parameters
+    integer, intent(in) :: nel, indt
+    real*8, intent(in) :: psi(indt, nel, 2), ainv(nel), ainvn(nel)
+    real*8, intent(inout) :: winv(nel, indt)
+
+    ! local variables
+    integer :: i, j
+
 #ifdef _OFFLOAD
     if (yes_ontarget) then
 !$omp target teams distribute parallel do collapse(2)
@@ -42,12 +50,48 @@ subroutine upwinvp(nel, indt, winv, ainv, ainvn, psi)
 #endif
     return
 end
+
+!> subroutine to update the matrix winv, for complex(16)
+subroutine upwinvp_complex(nel, indt, winv, ainv, ainvn, psi)
+    use constants, only: yes_ontarget
+    implicit none
+
+    ! argument parameters
+    integer, intent(in) :: nel, indt
+    complex*16, intent(in) :: psi(indt, nel, 2), ainv(nel), ainvn(nel)
+    complex*16, intent(inout) :: winv(nel, indt)
+
+    ! local variables
+    integer i, j
+
+#ifdef _OFFLOAD
+!$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
+#endif
+    do i = 1, indt
+        do j = 1, nel
+            winv(j, i) = winv(j, i) + ainv(j)*psi(i, j, 1) + ainvn(j)*psi(i, j, 2)
+        end do
+    end do
+#ifdef _OFFLOAD
+!$omp end target teams distribute parallel do
+#endif
+    return
+end
+
+!> subroutine to update the matrix winv of Pfaffian, for real(8)
 subroutine upwinvp_pfaff(nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt, winvup, winvdo, psi, ainv)
     use constants, only: yes_ontarget
     implicit none
-    integer nelc, nelcdo, nelup, neldo, nmol, nmolipf, nmolshift, indt, i, j
-    real*8 winvup(nelup, indt), winvdo(neldo, indt), psi(nmolipf, indt), ainv(nmol)
+
+    ! argument parameters
+    integer, intent(in) :: nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt
+    real*8, intent(in) :: psi(nmolipf, indt), ainv(nmol)
+    real*8, intent(inout) :: winvup(nelup, indt), winvdo(neldo, indt)
+
+    ! local variables
+    integer i, j, nelcdo
     real*8 csum
+
     if (yes_ontarget) then
     if (nelc .le. nelup) then
 #ifdef _OFFLOAD
@@ -116,11 +160,19 @@ subroutine upwinvp_pfaff(nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt, win
     end if
     return
 end
+
+!> subroutine to update the matrix winv of Pfaffian, for complex(16)
 subroutine upwinvp_pfaff_complex(nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt, winvup, winvdo, psi, ainv)
     use constants, only: yes_ontarget
     implicit none
-    integer nelc, nelcdo, nelup, neldo, nmol, nmolipf, nmolshift, indt, i, j
-    complex*16 winvup(nelup, indt), winvdo(neldo, indt), psi(nmolipf, indt), ainv(nmol)
+
+    ! argument parameters
+    integer, intent(in) :: nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt
+    complex*16, intent(in) :: psi(nmolipf, indt), ainv(nmol)
+    complex*16, intent(inout) :: winvup(nelup, indt), winvdo(neldo, indt)
+
+    ! local variables
+    integer i, j, nelcdo
     complex*16 csum
     if (yes_ontarget) then
         if (nelc .le. nelup) then
@@ -194,21 +246,3 @@ subroutine upwinvp_pfaff_complex(nelc, nelup, neldo, nmol, nmolipf, nmolshift, i
     end if
     return
 end
-subroutine upwinvp_complex(nel, indt, winv, ainv, ainvn, psi)
-    use constants, only: yes_ontarget
-    implicit none
-    integer nel, indt, i, j
-    complex*16 winv(nel, indt), psi(indt, nel, 2), ainv(nel), ainvn(nel)
-#ifdef _OFFLOAD
-!$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
-#endif
-    do i = 1, indt
-        do j = 1, nel
-            winv(j, i) = winv(j, i) + ainv(j)*psi(i, j, 1) + ainvn(j)*psi(i, j, 2)
-        end do
-    end do
-#ifdef _OFFLOAD
-!$omp end target teams distribute parallel do
-#endif
-    return
-end
diff --git a/src/m_common/write_type_orb.f90 b/src/m_common/write_type_orb.f90
index 1702557..0fbd79f 100644
--- a/src/m_common/write_type_orb.f90
+++ b/src/m_common/write_type_orb.f90
@@ -13,9 +13,15 @@
 ! You should have received a copy of the GNU General Public License
 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
 
+!> This subroutine returns typeorb for jastrow orbitals
 subroutine write_type_orb(nshellj, multij, ioccj, typeorb)
     implicit none
-    integer nshellj, multij(*), ioccj(*), typeorb(*)
+
+    ! argument parameters
+    integer, intent(in) :: nshellj, multij(*), ioccj(*)
+    integer, intent(out) :: typeorb(*)
+
+    ! local variables
     integer i, ii, j, ind_type
 
     ii = 0
diff --git a/src/m_common/zgemm_my.f90 b/src/m_common/zgemm_my.f90
index ba3baa0..571c505 100644
--- a/src/m_common/zgemm_my.f90
+++ b/src/m_common/zgemm_my.f90
@@ -13,6 +13,7 @@
 ! You should have received a copy of the GNU General Public License
 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
 
+!> This subroutine is a wrapper for the LAPACK zgemm routine
 subroutine ZGEMM_MY(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA&
         &, B, LDB, BETA, C, LDC, nproc, rank, comm_mpi)
     implicit none
diff --git a/src/m_common/zsktri.f90 b/src/m_common/zsktri.f90
index 4281420..855b12d 100644
--- a/src/m_common/zsktri.f90
+++ b/src/m_common/zsktri.f90
@@ -13,19 +13,29 @@
 ! You should have received a copy of the GNU General Public License
 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
 
+!> This subroutine manipulates a skewsymmetric matrix A.
 subroutine zsktri(uplo, n, a, lda, ainv, ldinv, ipiv, work, info)
     implicit none
-    character uplo
-    integer n, i, j, lda, ldinv, info, lwork, ipiv(*)
-    complex*16 a(lda, *), ainv(ldinv, *), work(*)
+
+    ! argument parameters
+    character, intent(in) :: uplo
+    integer, intent(in) :: n, lda, ldinv, ipiv(*)
+    complex*16, intent(in out) :: a(lda, *), work(*)
+    complex*16, intent(out) :: ainv(ldinv, *)
+    integer, intent(out) :: info
+
+    ! local variables
+    integer i, j
     real*8 rcond
     logical yeslap
+
     !     First factorization assumed
-    !      lwork=3*n
+    !     lwork=3*n
     !     CALL ZSKTRF( UPLO, 'N', N, A, LDA, IPIV, WORK, LWORK, INFO)
     !     Identity matrix
     !     ipiv dimension required 3n
     !     work dimension required n^2+12*n-2
+    info = 0
     yeslap = .false. ! if false the homemade algorithm is done.
     do i = 1, n
         do j = 1, i - 1
@@ -36,13 +46,13 @@ subroutine zsktri(uplo, n, a, lda, ainv, ldinv, ipiv, work, info)
             ainv(j, i) = dcmplx(0.d0, 0.d0)
         end do
     end do
-    !     fisrt permutation of ainv
+    ! fisrt permutation of ainv
     do i = n, 1, -1
         work(1:n) = ainv(i, 1:n)
         ainv(i, 1:n) = ainv(ipiv(i), 1:n)
         ainv(ipiv(i), 1:n) = work(1:n)
     end do
-    !      SECOND
+    ! second
     if (UPLO .eq. 'u' .or. UPLO .eq. 'U') then
         do i = 1, N - 1
             work(i) = a(i, i + 1)
@@ -63,7 +73,7 @@ subroutine zsktri(uplo, n, a, lda, ainv, ldinv, ipiv, work, info)
         a(2:n, 1) = dcmplx(0.d0, 0.d0)
     end if
     work(2*N - 1:3*N - 2) = dcmplx(0.d0, 0.d0) ! diagonal elements of skew matrix , obviously set to zero.
-    call ZTRTRS(UPLO, 'N', 'U', N, N, A, LDA, ainv, LDINV, INFO)
+    call ZTRTRS(UPLO, 'N', 'U', N, N, A, LDA, ainv, LDINV, INFO) ! ainv = A^-1
     ! USE STANDARD LAPACK DGTSVX or the simples DGTSV?
     if (yeslap) then
         call ZGTSVX('N', 'N', N, N, work(N), work(2*N - 1), work, work(3*N)&
@@ -83,7 +93,6 @@ subroutine zsktri(uplo, n, a, lda, ainv, ldinv, ipiv, work, info)
         ainv(i, 1:n) = ainv(ipiv(i), 1:n)
         ainv(ipiv(i), 1:n) = work(1:n)
     end do
-
     return
 end
 
diff --git a/src/m_common/zsktrs.f90 b/src/m_common/zsktrs.f90
index 2effe03..638a85e 100644
--- a/src/m_common/zsktrs.f90
+++ b/src/m_common/zsktrs.f90
@@ -1,4 +1,4 @@
-! Copyright (C) 2022 TurboRVB group
+! Copyright (C) 2023- TurboRVB group
 !
 ! This program is free software: you can redistribute it and/or modify
 ! it under the terms of the GNU General Public License as published by
@@ -13,11 +13,22 @@
 ! You should have received a copy of the GNU General Public License
 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
 
+!> This subroutine solves A X = B
+!> To be deleted!
 subroutine zsktrs(uplo, n, nhrs, a, b, ldb, x, info)
     implicit none
-    integer n, nhrs, ldb, info, i, j
-    complex*16 a(*), b(ldb, nhrs), x(nhrs, *)
-    character uplo
+
+    ! argument parameters
+    integer, intent(in) :: n, nhrs, ldb
+    integer, intent(out) :: info
+    complex*16, intent(in out) :: a(*), b(ldb, nhrs)
+    complex*16, intent(out) :: x(nhrs, *)
+    character, intent(in) :: uplo
+
+    ! local variables
+    integer i, j
+
+    info = 0
     if (uplo .eq. 'l' .or. uplo .eq. 'L') then
         a(1:n - 1) = -a(1:n - 1)
     end if
diff --git a/src/m_common/zsktrs_qp.f90 b/src/m_common/zsktrs_qp.f90
index aa0c9f9..272d47a 100644
--- a/src/m_common/zsktrs_qp.f90
+++ b/src/m_common/zsktrs_qp.f90
@@ -13,16 +13,31 @@
 ! You should have received a copy of the GNU General Public License
 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
 
+!> This subroutine solves A * X = B.
+!> Our homemade algorithm to solve A * X = B (No transpose)
+!> where A is a skewsymmetric matrix, X and B are rectanglar matrices.
+!> The input vector a stores the elements of A such that
+!> A(i+1,i)=a(i) and A(i,i+1)=-a(i) for i=1,..,n-1
+!> a(n), b(ldb, n), x(n)
 subroutine zsktrs(uplo, n, nhrs, a, b, ldb, x, info)
     implicit none
-    integer n, nhrs, ldb, info, i, j
-    complex*16 a(*), b(ldb, *), x(*)
+
+    ! argument variables
+    integer, intent(in) :: n, nhrs, ldb
+    integer, intent(out) :: info
+    complex*16, intent(inout) :: a(*), b(ldb, *)
+    complex*16, intent(out) :: x(*)
+    character, intent(in) :: uplo
+
+    ! local variables
+    integer i, j
 #ifdef __PORT
     complex*16 aq, a2q, xo
 #else
     complex*32 aq, a2q, xo
 #endif
-    character uplo
+
+    info = 0
     if (uplo .eq. 'l' .or. uplo .eq. 'L') then
         a(1:n - 1) = -a(1:n - 1)
     end if