diff --git a/cmake/Sources.cmake b/cmake/Sources.cmake
index dea3d26a..e5ec00c4 100644
--- a/cmake/Sources.cmake
+++ b/cmake/Sources.cmake
@@ -80,6 +80,7 @@ set(IRON_Fortran_SRC
     generated_mesh_routines.F90
     generated_mesh_access_routines.F90
     Hamilton_Jacobi_equations_routines.F90
+    hash_routines.F90
     Helmholtz_equations_routines.F90
     #Helmholtz_TEMPLATE_equations_routines.F90
     history_routines.F90
diff --git a/src/basis_routines.F90 b/src/basis_routines.F90
index 98f14871..b26d709e 100644
--- a/src/basis_routines.F90
+++ b/src/basis_routines.F90
@@ -6936,7 +6936,7 @@ SUBROUTINE Gauss_Simplex(order,numberOfVertices,n,x,w,err,error,*)
           CALL FlagError(localError,err,error,*999)
         ENDIF
         !Gauss point 1
-        x(1,1)=1.0_DP/2.0
+        x(1,1)=1.0_DP/2.0_DP
         x(2,1)=1.0_DP/2.0_DP
         w(1)=1.0_DP
       CASE(2)
@@ -6993,7 +6993,7 @@ SUBROUTINE Gauss_Simplex(order,numberOfVertices,n,x,w,err,error,*)
         ENDIF
         !Gauss point 1
         x(1,1)=(1.0_DP-SQRT(0.6_DP))/2.0_DP
-        x(2,1)=1-x(1,1)
+        x(2,1)=1.0_DP-x(1,1)
         w(1)=5.0_DP/18.0_DP
         !Gauss point 2
         x(1,2)=1.0_DP/2.0_DP
@@ -7017,7 +7017,7 @@ SUBROUTINE Gauss_Simplex(order,numberOfVertices,n,x,w,err,error,*)
         ENDIF
         !Gauss point 1
         x(1,1)=(1.0_DP-SQRT(0.6_DP))/2.0_DP
-        x(2,1)=1-x(1,1)
+        x(2,1)=1.0_DP-x(1,1)
         w(1)=5.0_DP/18.0_DP
         !Gauss point 2
         x(1,2)=1.0_DP/2.0_DP
@@ -7099,8 +7099,8 @@ SUBROUTINE Gauss_Simplex(order,numberOfVertices,n,x,w,err,error,*)
           CALL FlagError(localError,err,error,*999)
         ENDIF
         lC=1.0_DP/3.0_DP
-        wC=-9.0_DP/16.0
-        alpha1=25.0_DP
+        wC=-3.0_DP/4.0_DP
+        alpha1=2.0_DP/5.0_DP
         wAlpha1=25.0_DP/48.0_DP
         l1Alpha1=(1.0_DP+2.0_DP*alpha1)/3.0_DP
         l2Alpha1=(1.0_DP-alpha1)/3.0_DP
diff --git a/src/boundary_condition_routines.F90 b/src/boundary_condition_routines.F90
index 2a9e9368..4c2562aa 100644
--- a/src/boundary_condition_routines.F90
+++ b/src/boundary_condition_routines.F90
@@ -325,8 +325,16 @@ SUBROUTINE BOUNDARY_CONDITIONS_CREATE_FINISH(BOUNDARY_CONDITIONS,ERR,ERROR,*)
                   IF(ASSOCIATED(VARIABLE_DOMAIN_MAPPING)) THEN
                     SEND_COUNT=VARIABLE_DOMAIN_MAPPING%NUMBER_OF_GLOBAL
                     IF(computationalEnvironment%numberOfComputationalNodes>1) THEN
-                      !\todo This operation is a little expensive as we are doing an unnecessary sum across all the ranks in order to combin
+                      !\todo This operation is a little expensive as we are doing an unnecessary sum across all the ranks in order to combine
                       !\todo the data from each rank into all ranks. We will see how this goes for now.
+
+                      ! Should be replaced with scatter(?)!! And local implementation of dof_types!!!
+                      ! But solver mapping STILL requires dof types global
+                      ! Can change this but would work only here.
+                      ! ALLREDUCE must stay until solver_mapping is not changed.
+               
+
+                      ! Summing up all global vectors returns the number of fixed dofs in the row (?)
                       CALL MPI_ALLREDUCE(MPI_IN_PLACE,BOUNDARY_CONDITION_VARIABLE%DOF_TYPES, &
                         & SEND_COUNT,MPI_INTEGER,MPI_SUM,computationalEnvironment%mpiCommunicator,MPI_IERROR)
                       CALL MPI_ERROR_CHECK("MPI_ALLREDUCE",MPI_IERROR,ERR,ERROR,*999)
@@ -408,8 +416,12 @@ SUBROUTINE BOUNDARY_CONDITIONS_CREATE_FINISH(BOUNDARY_CONDITIONS,ERR,ERROR,*)
                     BOUNDARY_CONDITIONS_DIRICHLET=>BOUNDARY_CONDITION_VARIABLE%DIRICHLET_BOUNDARY_CONDITIONS
                     IF(ASSOCIATED(BOUNDARY_CONDITIONS_DIRICHLET)) THEN
                       ! Find dirichlet conditions
+
+                      ! replace with (total) n. of local!!??
                       dirichlet_idx=1
                       DO dof_idx=1,FIELD_VARIABLE%NUMBER_OF_GLOBAL_DOFS
+                        
+                        ! BC_DOF_FIXED is 1
                         IF(BOUNDARY_CONDITION_VARIABLE%DOF_TYPES(dof_idx)==BOUNDARY_CONDITION_DOF_FIXED) THEN
                           BOUNDARY_CONDITIONS_DIRICHLET%DIRICHLET_DOF_INDICES(dirichlet_idx)=dof_idx
                           dirichlet_idx=dirichlet_idx+1
@@ -1705,6 +1717,7 @@ SUBROUTINE BoundaryConditions_SetConditionType(boundaryConditionsVariable,global
       END IF
     END IF
     !Update Dirichlet DOF count
+    ! should be local!!!
     previousDof=boundaryConditionsVariable%DOF_TYPES(globalDof)
     IF(dofType==BOUNDARY_CONDITION_DOF_FIXED.AND.previousDof/=BOUNDARY_CONDITION_DOF_FIXED) THEN
       boundaryConditionsVariable%NUMBER_OF_DIRICHLET_CONDITIONS= &
@@ -1716,6 +1729,7 @@ SUBROUTINE BoundaryConditions_SetConditionType(boundaryConditionsVariable,global
 
     !Set the boundary condition type and DOF type
     boundaryConditionsVariable%CONDITION_TYPES(globalDof)=condition
+    ! should be local!!!
     boundaryConditionsVariable%DOF_TYPES(globalDof)=dofType
     IF(DIAGNOSTICS1) THEN
       CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"Boundary Condition Being Set",err,error,*999)
@@ -2325,6 +2339,7 @@ SUBROUTINE BoundaryConditions_NeumannMatricesInitialise(boundaryConditionsVariab
 
         CALL DOMAIN_MAPPINGS_LOCAL_FROM_GLOBAL_CALCULATE(pointDofMapping,err,error,*999)
 
+        ! here allocation of integrationMatrix!!!???
         CALL DistributedMatrix_CreateStart(rowMapping,pointDofMapping,boundaryConditionsNeumann%integrationMatrix,err,error,*999)
         SELECT CASE(boundaryConditionsVariable%BOUNDARY_CONDITIONS%neumannMatrixSparsity)
         CASE(BOUNDARY_CONDITION_SPARSE_MATRICES)
@@ -2500,6 +2515,8 @@ SUBROUTINE BoundaryConditions_NeumannMatricesInitialise(boundaryConditionsVariab
 
           CALL DistributedMatrix_StorageTypeSet(boundaryConditionsNeumann%integrationMatrix, &
             & DISTRIBUTED_MATRIX_COMPRESSED_ROW_STORAGE_TYPE,err,error,*999)
+
+          ! replace with local quantities!!!???
           CALL DistributedMatrix_NumberOfNonZerosSet(boundaryConditionsNeumann%integrationMatrix,numberNonZeros,err,error,*999)
           CALL DistributedMatrix_StorageLocationsSet(boundaryConditionsNeumann%integrationMatrix, &
             & rowIndices,columnIndices(1:numberNonZeros),err,error,*999)
@@ -2724,6 +2741,7 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,*
             neumannNodeNumber=rhsVariable%DOF_TO_PARAM_MAP%NODE_DOF2PARAM_MAP(3,neumannDofNyy)
             SELECT CASE(rhsVariable%COMPONENTS(componentNumber)%DOMAIN%NUMBER_OF_DIMENSIONS)
             CASE(1)
+              ! replace with local????
               CALL DistributedMatrix_ValuesSet(neumannConditions%integrationMatrix,neumannLocalDof,neumannDofIdx, &
                 & 1.0_DP,err,error,*999)
             CASE(2)
@@ -2816,6 +2834,7 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,*
                     END IF
 
                     ! Add integral term to N matrix
+                    ! Already local, ok???
                     CALL DistributedMatrix_ValuesAdd(neumannConditions%integrationMatrix,localDof,neumannDofIdx, &
                       & integratedValue,err,error,*999)
                   END DO
@@ -2911,6 +2930,7 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,*
                     END IF
 
                     ! Add integral term to N matrix
+                    ! Already local, ok???
                     CALL DistributedMatrix_ValuesAdd(neumannConditions%integrationMatrix,localDof,neumannDofIdx, &
                       & integratedValue,err,error,*999)
                   END DO
@@ -2944,6 +2964,8 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,*
         & integratedValues,err,error,*999)
       CALL DistributedVector_AllValuesSet(integratedValues,0.0_DP,err,error,*999)
       ! Perform matrix multiplication, f = N q, to calculate force vector from integration matrix and point values
+
+      ! local??? Mapping on f???
       CALL DistributedMatrix_MatrixByVectorAdd(DISTRIBUTED_MATRIX_VECTOR_NO_GHOSTS_TYPE,1.0_DP, &
         & neumannConditions%integrationMatrix,neumannConditions%pointValues,integratedValues, &
         & err,error,*999)
@@ -3251,6 +3273,7 @@ SUBROUTINE BoundaryConditions_DofConstraintSet(boundaryConditions,fieldVariable,
     END DO
 
     !Check DOFs are free
+    ! Should be local!
     DO dofIdx=1,numberOfDofs
       IF(boundaryConditionsVariable%dof_types(dofs(dofIdx))/=BOUNDARY_CONDITION_DOF_FREE) THEN
         CALL FlagError("DOF number "//TRIM(NumberToVstring(dofs(dofIdx),"*",err,error))// &
@@ -3375,6 +3398,8 @@ SUBROUTINE BoundaryConditions_DofConstraintsCreateFinish(boundaryConditionsVaria
           !Check that the constrained DOFs are still set to be constrained, as
           !subsequently setting a boundary condition would change the DOF type but
           !not update the DOF constraints structure.
+
+          ! should be local!!
           IF(boundaryConditionsVariable%dof_types(globalDof)/=BOUNDARY_CONDITION_DOF_CONSTRAINED) THEN
             CALL FlagError("Global DOF number "//TRIM(NumberToVstring(globalDof,"*",err,error))// &
               & " is part of a linear constraint but the DOF type has been changed"// &
@@ -3385,6 +3410,7 @@ SUBROUTINE BoundaryConditions_DofConstraintsCreateFinish(boundaryConditionsVaria
             globalDof2=dofConstraint%dofs(dofIdx)
             localDof2=variableDomainMapping%global_to_local_map(globalDof2)%local_number(1)
             !Check a Dirichlet conditions hasn't also been set on this DOF
+            ! should be local!!
             IF(boundaryConditionsVariable%dof_types(globalDof2)/=BOUNDARY_CONDITION_DOF_FREE) THEN
               CALL FlagError("A Dirichlet boundary condition has been set on DOF number "// &
                 & TRIM(NumberToVstring(globalDof2,"*",err,error))// &
@@ -3684,10 +3710,15 @@ SUBROUTINE BOUNDARY_CONDITIONS_VARIABLE_INITIALISE(BOUNDARY_CONDITIONS,FIELD_VAR
             BOUNDARY_CONDITIONS_VARIABLE%VARIABLE=>FIELD_VARIABLE
             ALLOCATE(BOUNDARY_CONDITIONS_VARIABLE%CONDITION_TYPES(VARIABLE_DOMAIN_MAPPING%NUMBER_OF_GLOBAL),STAT=ERR)
             IF(ERR/=0) CALL FlagError("Could not allocate global boundary condition types.",ERR,ERROR,*999)
+
+            ! THIS ALLOCATION SHOULD BE LOCAL!
             ALLOCATE(BOUNDARY_CONDITIONS_VARIABLE%DOF_TYPES(VARIABLE_DOMAIN_MAPPING%NUMBER_OF_GLOBAL),STAT=ERR)
             IF(ERR/=0) CALL FlagError("Could not allocate global boundary condition dof types.",ERR,ERROR,*999)
             BOUNDARY_CONDITIONS_VARIABLE%CONDITION_TYPES=BOUNDARY_CONDITION_FREE
+            
+            ! All initialised to zero!!!
             BOUNDARY_CONDITIONS_VARIABLE%DOF_TYPES=BOUNDARY_CONDITION_DOF_FREE
+            
             ALLOCATE(BOUNDARY_CONDITIONS_VARIABLE%DOF_COUNTS(MAX_BOUNDARY_CONDITION_NUMBER),STAT=ERR)
             IF(ERR/=0) CALL FlagError("Could not allocate boundary condition DOF counts array.",ERR,ERROR,*999)
             BOUNDARY_CONDITIONS_VARIABLE%DOF_COUNTS=0
diff --git a/src/distributed_matrix_vector.F90 b/src/distributed_matrix_vector.F90
index 927e4a10..7dcf4ee2 100644
--- a/src/distributed_matrix_vector.F90
+++ b/src/distributed_matrix_vector.F90
@@ -60,6 +60,7 @@ MODULE DistributedMatrixVector
 #endif
   USE Strings
   USE Types
+  USE HashRoutines
   USE LINKEDLIST_ROUTINES
 
 #include "macros.h"
@@ -974,7 +975,22 @@ SUBROUTINE DistributedMatrix_CMISSInitialise(distributedMatrix,err,error,*)
     INTEGER(INTG) :: dummyErr
     TYPE(DOMAIN_MAPPING_TYPE), POINTER :: rowDomainMapping,columnDomainMapping
     TYPE(VARYING_STRING) :: dummyError,localError
-    
+
+    ! Hash variables
+    INTEGER(INTG), ALLOCATABLE :: SKey(:) 
+    INTEGER(INTG), ALLOCATABLE :: SIntVal(:,:)
+    REAL(SP), ALLOCATABLE      :: SRealVal(:,:)
+
+    INTEGER(INTG) :: i, n, q, myComputationalNodeNumber, indexFound, dataType, dataSize, newDataSize
+    LOGICAL :: isFound
+
+    INTEGER(INTG), ALLOCATABLE :: valueInt(:) ! Allow for different returned value INT sizes
+    REAL(SP), ALLOCATABLE      :: valueSp(:)  ! Allow for different returned value SP sizes
+
+    INTEGER(INTG), ALLOCATABLE :: valueIntMult(:,:) ! Same in case of multiple lists
+    REAL(SP),      ALLOCATABLE :: valueSpMult (:,:) ! Same in case of multiple lists
+    INTEGER(INTG), ALLOCATABLE :: listNum(:) ! Array of list num for getting values from many lists 
+
     ENTERS("DistributedMatrix_CMISSInitialise",err,error,*998)
 
     IF(.NOT.ASSOCIATED(distributedMatrix)) CALL FlagError("Distributed matrix is not associated.",err,error,*998)
@@ -1009,7 +1025,171 @@ SUBROUTINE DistributedMatrix_CMISSInitialise(distributedMatrix,err,error,*)
         & TRIM(NumberToVString(distributedMatrix%ghostingType,"*",err,error))//" is invalid."
       CALL FlagError(localError,err,error,*999)
     END SELECT
-    
+
+    ! START my hash implementation    
+
+    ! use columnDomainMapping ltg info to create hash:
+    ! hash type in type of distributedMatrix%cmiss, move in future? petsc already handles. OK
+
+    ! Decide data type of values to add to the table
+    dataType = LIST_SP_TYPE 
+    ! dataType = LIST_INTG_TYPE
+    ! Decide data size of values
+    dataSize = 3
+
+    myComputationalNodeNumber = ComputationalEnvironment_NodeNumberGet(ERR,ERROR)
+
+    IF (diagnostics1) THEN
+      CALL WriteStringVector(DIAGNOSTIC_OUTPUT_TYPE, 1,1, size(columnDomainMapping%LOCAL_TO_GLOBAL_MAP,1), &
+      & 4,4, columnDomainMapping%LOCAL_TO_GLOBAL_MAP, & 
+      & '("Printing out the map  :", 4(x,I8))','(23x,4(x,I8))', ERR, ERROR, *999)
+    END IF
+
+    n = size(columnDomainMapping%LOCAL_TO_GLOBAL_MAP)
+
+    ! Create a new table
+    ! If a table already exists, content is zeroed at CreateStart. 
+    NULLIFY(distributedMatrix%cmiss%columnHashTable)
+
+    CALL HashTable_CreateStart(distributedMatrix%cmiss%columnHashTable,ERR,ERROR,*999)
+    ! define some parameters here if needed
+    CALL HashTable_CreateFinish(distributedMatrix%cmiss%columnHashTable,ERR,ERROR,*999)
+
+    ! Create array of keys
+    ALLOCATE(SKey(n),STAT=err)
+    IF(ERR/=0) CALL FlagError("Could not allocate array",ERR,ERROR,*999)
+    SKey = columnDomainMapping%LOCAL_TO_GLOBAL_MAP(1:n)
+
+    ! Create initial arrays of values
+    ALLOCATE(SIntVal(dataSize,n),STAT=err)
+    IF(ERR/=0) CALL FlagError("Could not allocate array",ERR,ERROR,*999)
+    DO i=1,dataSize  ! Just repeat the same (nonsense) info (= local numbering) 
+      SIntVal(i,1:n)=[ (i, i=1,n)] +100
+    END DO
+
+    ALLOCATE(SRealVal(dataSize,n),STAT=err)
+    IF(ERR/=0) CALL FlagError("Could not allocate array",ERR,ERROR,*999)
+    DO i=1,dataSize  ! Just repeat the same (nonsense) info (= local numbering + real part)
+      SRealVal(i,1:n) = [ (i, i=1,n)] +100.15_SP
+    END DO
+
+    SELECT CASE (dataType)
+
+    CASE(LIST_INTG_TYPE)
+
+    ! Create the list of values and compute the table based on the key array
+    ! Create a new table
+    ! .FALSE. = new table (no addition)
+    CALL HashTable_ValuesSetAndInsert(distributedMatrix%cmiss%columnHashTable, SKey(1:n-3), &
+      & SIntVal(1:dataSize,1:n-3), .FALSE., ERR, ERROR, *999)
+
+    ! Test insertion to already created table
+    ! Recompute the table adding new keys and expand the list of values
+    ! .TRUE. = add new elements
+    CALL HashTable_ValuesSetAndInsert(distributedMatrix%cmiss%columnHashTable, SKey(n-2:n), &
+      & SIntVal(1:dataSize,n-2:n), .TRUE., ERR, ERROR, *999)
+
+    CASE(LIST_SP_TYPE)
+
+    ! Create the list of values and compute the table based on the key array
+    ! Create a new table
+    ! .FALSE. = new table (no addition)
+    CALL HashTable_ValuesSetAndInsert(distributedMatrix%cmiss%columnHashTable, SKey(1:n-3), &
+      & SRealVal(1:dataSize,1:n-3), .FALSE., ERR, ERROR, *999)
+
+    ! Test insertion to already created table
+    ! Recompute the table adding new keys and expand the list of values
+    ! .TRUE. = add new elements
+    CALL HashTable_ValuesSetAndInsert(distributedMatrix%cmiss%columnHashTable, SKey(n-2:n), &
+      & SRealVal(1:dataSize,n-2:n), .TRUE., ERR, ERROR, *999)
+
+    CASE DEFAULT
+      CALL FlagError("Invalid data type!!!",err,error,*999)  
+    END SELECT
+
+
+
+    ! Test query for q
+    IF (diagnostics1) THEN
+      DO q=1,50
+      ! Get the INDEX in the input vector (Skey, SIntVal) where q is found in the table (+ flag isFound)
+        CALL HashTable_GetKey(distributedMatrix%cmiss%columnHashTable, q, indexFound, isFound, ERR, ERROR, *999)
+        IF (isFound) THEN
+          SELECT CASE (dataType)
+          CASE(LIST_INTG_TYPE)
+            CALL HashTable_GetValue(distributedMatrix%cmiss%columnHashTable, indexFound, valueInt, ERR, ERROR, *999)
+ !           PRINT *, q, "Found! With value ", valueInt
+            CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "Key found: ", q, err, error, *999)
+            CALL WriteStringVector(DIAGNOSTIC_OUTPUT_TYPE, 1,1, size(valueInt,1), &
+              & 4,4, valueInt, & 
+              & '("with value     :", 4(x,I8))','(5x,4(x,I8))', ERR, ERROR, *999)
+          CASE(LIST_SP_TYPE)
+            CALL HashTable_GetValue(distributedMatrix%cmiss%columnHashTable, indexFound, valueSp, ERR, ERROR, *999)
+!            PRINT *, q, "Found! With value ", valueSp
+            CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "Key found: ", q, err, error, *999)
+            CALL WriteStringVector(DIAGNOSTIC_OUTPUT_TYPE, 1,1, size(valueSp,1), &
+              & 4,4, valueSp, & 
+              & '("with value     :", 4(x,F8.4))','(5x,4(x,F8.4))', ERR, ERROR, *999)
+          CASE DEFAULT
+            CALL FlagError("Invalid data type!!!",err,error,*999)  
+          END SELECT
+        ELSE
+          CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "Key does not exist! ", q, err, error, *999)
+        END IF 
+      END DO
+    END IF
+
+  ! Add additional values to the existing ones in table%ListSVal as
+  ! table%arrayOfListSVal(i)%ptr
+  ! Can be called AFTER the table has been created (table%ListSVal should already exist)
+  ! Add INTEGER or REAL values of different dimension
+    newDataSize = 2
+  ! List 2
+    CALL HashTable_AdditionalValuesSet (distributedMatrix%cmiss%columnHashTable, &
+      & SRealVal(1:newDataSize+1,1:n), LIST_SP_TYPE, err,error,*999)
+  ! List 3
+    CALL HashTable_AdditionalValuesSet (distributedMatrix%cmiss%columnHashTable, &
+      & SIntVal(1:newDataSize,1:n), LIST_INTG_TYPE, err,error,*999)
+  ! List 4
+    CALL HashTable_AdditionalValuesSet (distributedMatrix%cmiss%columnHashTable, &
+      & SRealVal(1:newDataSize,1:n), LIST_SP_TYPE, err,error,*999)
+
+
+! Test query in multiple values case for 1 key
+! The obtained value (real or intg) is sized as #lists X max dimension of the data in all lists
+! Then, there will be zeros if the value has smaller dimension. (Fix??)
+    IF (diagnostics1) THEN
+      q = 25
+      CALL HashTable_GetKey(distributedMatrix%cmiss%columnHashTable, q, indexFound, isFound, ERR, ERROR, *999)
+      IF (isFound) THEN
+        CALL HashTable_GetValue(distributedMatrix%cmiss%columnHashTable, indexFound, valueSpMult, listNum, ERR, ERROR, *999)
+        DO i=1,size(listNum,1)
+          !PRINT *, q, "Found Real!    With value ", valueSpMult(i,:), "in list number", listNum(i)
+
+          CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "Found in list number ", listNum(i), err, error, *999)
+          CALL WriteStringVector(DIAGNOSTIC_OUTPUT_TYPE, 1,1, size(valueSpMult(i,:),1), &
+            & 4,4, valueSpMult(i,:), & 
+            & '("REAL SP value     :", 4(x,F8.4))','(5x,4(x,F8.4))', ERR, ERROR, *999)
+
+        END DO
+        CALL HashTable_GetValue(distributedMatrix%cmiss%columnHashTable, indexFound, valueIntMult, listNum, ERR, ERROR, *999)
+        DO i=1,size(listNum,1)
+          !PRINT *, q, "Found Integer! With value ", valueIntMult(i,:), "in list number", listNum(i)
+          CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "Found in list number: ", listNum(i), err, error, *999)
+          CALL WriteStringVector(DIAGNOSTIC_OUTPUT_TYPE, 1,1, size(valueIntMult(i,:),1), &
+            & 4,4, valueIntMult(i,:), & 
+            & '("INTEGER value     :", 4(x,I8))','(5x,4(x,I8))', ERR, ERROR, *999)
+        END DO
+      
+      ELSE
+        CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "Key does not exist! ", q, err, error, *999)
+      END IF
+    END IF
+
+    STOP
+
+! end my hash_implementation
+
     EXITS("DistributedMatrix_CMISSInitialise")
     RETURN
 999 IF(ASSOCIATED(distributedMatrix%cmiss)) &
@@ -2663,7 +2843,8 @@ SUBROUTINE DistributedMatrix_PETScInitialise(distributedMatrix,err,error,*)
     distributedMatrix%petsc%useOverrideMatrix=.FALSE.
     CALL Petsc_MatInitialise(distributedMatrix%petsc%matrix,err,error,*999)
     CALL Petsc_MatInitialise(distributedMatrix%petsc%overrideMatrix,err,error,*999)
-   
+
+
     EXITS("DistributedMatrix_PETScInitialise")
     RETURN
 999 IF(ASSOCIATED(distributedMatrix%petsc)) &
diff --git a/src/hash_routines.F90 b/src/hash_routines.F90
new file mode 100644
index 00000000..2e0e390f
--- /dev/null
+++ b/src/hash_routines.F90
@@ -0,0 +1,1770 @@
+! Implements hash tables for static sets (i.e. do not change) as described in
+! Fredman, Komlos and Szemeredi 1984
+! using (arrays of) Lists that allow for INTG or REAL (implemented SP) of different sizes as data content (values)
+
+MODULE HashRoutines
+
+  USE BaseRoutines
+  USE INPUT_OUTPUT
+  USE CONSTANTS
+  USE KINDS
+  USE ISO_VARYING_STRING
+  USE STRINGS 
+  USE ComputationEnvironment
+  USE Types
+  USE Lists
+
+#include "macros.h"  
+
+  IMPLICIT NONE
+  
+  PRIVATE  ! by default
+
+  !Module parameters
+  INTEGER(INTG), PARAMETER :: HASH_INITIAL_ADDITIONAL_VALUES_SIZE=5 !<Initial number of additional lists of values 
+
+  ! TYPEs
+
+  ! Items for Linked Lists with 2 entries
+  TYPE HashLinkedListItemType
+    PRIVATE
+    INTEGER(INTG) :: key    
+    INTEGER(INTG) :: value
+    TYPE(HashLinkedListItemType),POINTER :: next => NULL()
+  END TYPE
+
+  ! Linked Lists
+  TYPE HashLinkedListType
+    PRIVATE
+    TYPE(HashLinkedListItemType),POINTER :: root => NULL()
+    TYPE(HashLinkedListItemType),POINTER :: last => NULL()
+  END TYPE
+
+  TYPE HashLinkedListPtrType
+    PRIVATE
+    TYPE(HashLinkedListType), POINTER :: ptr
+  END TYPE HashLinkedListPtrType
+
+  ! Array of Linked Lists
+  TYPE HashListArrayType
+     PRIVATE
+     TYPE(HashLinkedListPtrType), ALLOCATABLE :: vec(:)
+     INTEGER(INTG)                            :: vecLen
+     LOGICAL(INTG)                            :: isFinished
+  END TYPE HashListArrayType
+
+  ! The Hash Table
+  ! TYPE HashTableType
+  ! in types.f90
+
+  INTERFACE HashLinkedList_Add
+    MODULE PROCEDURE HashLinkedList_AddData
+    MODULE PROCEDURE HashLinkedList_AddList
+  END INTERFACE HashLinkedList_Add
+
+  INTERFACE HashTable_ListValuesSet
+    MODULE PROCEDURE HashTable_ListRealSpValuesSet
+    MODULE PROCEDURE HashTable_ListIntValuesSet
+  END INTERFACE HashTable_ListValuesSet
+
+  INTERFACE HashTable_ValuesSetAndInsert
+    MODULE PROCEDURE HashTable_IntValuesSetAndInsert
+    MODULE PROCEDURE HashTable_RealSpValuesSetAndInsert
+  END INTERFACE HashTable_ValuesSetAndInsert
+
+  INTERFACE HashTable_GetValue
+    MODULE PROCEDURE HashTable_GetRealSpValue
+    MODULE PROCEDURE HashTable_GetIntValue
+    MODULE PROCEDURE HashTable_GetListArrayIntValue
+    MODULE PROCEDURE HashTable_GetListArrayRealSpValue
+  END INTERFACE HashTable_GetValue
+
+  INTERFACE HashTable_AdditionalValuesSet
+    MODULE PROCEDURE HashTable_AdditionalIntValuesSet
+    MODULE PROCEDURE HashTable_AdditionalRealSpValuesSet
+  END INTERFACE HashTable_AdditionalValuesSet
+
+  ! public subs
+  ! Array of lists (keep private)
+  !PUBLIC :: HashListArray_CreateStart, HashListArray_CreateFinish
+  ! Table create
+  PUBLIC :: HashTable_ValuesSetAndInsert, HashTable_CreateStart, HashTable_CreateFinish, HashTable_AdditionalValuesSet
+  ! Table get values
+  PUBLIC :: HashTable_GetKey, HashTable_GetValue
+
+CONTAINS
+
+!==================================================
+
+  !> initialises or adds a piece of data to list
+  SUBROUTINE HashLinkedList_AddData(list,key, value, isDuplicateKey, err,error,*)
+
+    !Argument Variables
+    TYPE(HashLinkedListType), POINTER :: list
+
+    INTEGER(INTG),INTENT(IN) :: key, value
+
+    INTEGER(INTG), INTENT(OUT) :: err !<The error code
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    LOGICAL, INTENT (OUT) :: isDuplicateKey
+
+    !Local Variables
+    TYPE(HashLinkedListItemType),POINTER :: current => NULL()
+
+    ENTERS("HashLinkedList_AddData",err,error,*999)
+
+    isDuplicateKey = .FALSE.
+
+    ! Allocate the list if not done already
+    IF(.NOT. ASSOCIATED(list)) THEN
+        ALLOCATE(list, STAT=err)
+        IF(err/=0) CALL FlagError("Could not allocate list",err,error,*999)
+    END IF    
+
+    IF(ASSOCIATED(list%root)) THEN
+      current => list%root
+       DO
+         ! key is there already
+         IF(current%key==key) THEN
+         isDuplicateKey = .TRUE.
+         EXIT ! do nothing
+         END IF 
+
+         ! tail is reached
+         IF(.NOT. ASSOCIATED(current%next)) THEN
+            ALLOCATE(current%next)
+            current%next%key   = key
+            current%next%value = value
+            list%last => current%next
+         EXIT
+         END IF
+          
+         ! proceed
+         current => current%next 
+
+        END DO  
+       
+     ELSE     ! list is empty: create root 
+
+      ALLOCATE(list%root)
+      list%root%key   = key
+      list%root%value = value
+      list%last => list%root
+    END IF
+
+    EXITS("HashLinkedList_AddData")
+    RETURN
+999 ERRORSEXITS("HashLinkedList_AddData",err,error)
+    RETURN 1
+
+  END SUBROUTINE HashLinkedList_AddData
+
+!==================================================
+
+  !> Returns length of list
+  SUBROUTINE HashLinkedList_Size(list,n,err,error,*)
+
+    ! argument variables
+    TYPE(HashLinkedListType), POINTER :: list
+    INTEGER(INTG), INTENT(OUT) :: n
+    INTEGER(INTG), INTENT(OUT) :: err !<The error code
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+ 
+    ! local variables
+    TYPE(HashLinkedListItemType),POINTER :: current => NULL()
+    LOGICAL :: isEmpty
+
+    ENTERS("HashLinkedList_Size", err, error, *999)
+
+    ! Check if list is allocated
+    ! Return zero-size array also if list is empty
+    CALL HashLinkedList_IsEmpty(list, isEmpty, err, error, *999)
+    IF((.NOT. ASSOCIATED(list)).OR. isEmpty) THEN
+        n = 0
+    ELSE
+
+    ! Traversing to find size
+     current => list%root
+     n=1
+     DO
+       IF(ASSOCIATED(current%next)) THEN
+         n=n+1
+         current => current%next
+       ELSE
+         EXIT
+       END IF
+     END DO
+    
+    END IF
+
+
+    EXITS("HashLinkedList_Size")
+    RETURN
+999 ERRORSEXITS("HashLinkedList_Size",err,error)
+    RETURN 1
+
+  END SUBROUTINE HashLinkedList_Size
+
+!==================================================
+
+  SUBROUTINE HashLinkedList_GetData(list,key, value, isFound, err, error, *)
+
+    ! argument variables
+    TYPE(HashLinkedListType), POINTER :: list
+    INTEGER(INTG),INTENT(IN)   :: key
+    INTEGER(INTG),INTENT(OUT)  :: value
+    INTEGER(INTG), INTENT(OUT) :: err !<The error code
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+
+    ! local variables
+    TYPE(HashLinkedListItemType),POINTER :: current => NULL()
+    LOGICAL :: isFound
+
+    ENTERS("HashLinkedList_GetData", err, error, *999)
+ 
+    isFound = .FALSE.
+   
+    IF(ASSOCIATED(list%root)) THEN
+    
+    current => list%root
+    DO
+         ! key is found 
+         IF(current%key==key)THEN
+          value = current%value 
+           
+          isFound = .TRUE.
+
+         END IF 
+
+         ! tail is reached
+         IF(.NOT. ASSOCIATED(current%next)) EXIT
+          
+         ! otherwise proceed
+         current => current%next 
+    END DO
+  
+    END IF
+   
+    EXITS("HashLinkedList_GetData")
+    RETURN
+999 ERRORSEXITS("HashLinkedList_GetData",err,error)
+    RETURN 1
+
+
+   END SUBROUTINE HashLinkedList_GetData
+
+!==================================================
+
+  !> adds all data from one list (addList) to another (list)
+  SUBROUTINE HashLinkedList_AddList(list,addList,err,error, *)
+
+    ! argument variables
+    TYPE(HashLinkedListType), POINTER :: list
+    TYPE(HashLinkedListType), POINTER :: addList
+    INTEGER(INTG), INTENT(OUT) :: err !<The error code
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    ! local variables
+    LOGICAL :: isFound, isEmpty
+    TYPE(HashLinkedListItemType),POINTER :: current => NULL()
+
+    ENTERS("HashLinkedList_AddList", err, error, *999)
+
+    CALL HashLinkedList_IsEmpty(addList, isEmpty, err, error, *999)
+    IF(.NOT. isEmpty) THEN
+
+    current => addList%root
+    DO
+      CALL HashLinkedList_AddData(list,current%key, current%value, isFound, err,error,*999)
+      IF(ASSOCIATED(current%next)) THEN
+        current => current%next
+      ELSE
+        EXIT
+      END IF
+    END DO
+
+    END IF
+
+    EXITS("HashLinkedList_AddList")
+    RETURN
+999 ERRORSEXITS("HashLinkedList_AddList",err,error)
+    RETURN 1
+
+  END SUBROUTINE HashLinkedList_AddList
+
+!==================================================
+
+  !> removes the first item from list and returns its value in data
+  SUBROUTINE HashLinkedList_RemoveFirst(list,key,value,err,error,*)
+ 
+    ! Argument Variables
+    TYPE(HashLinkedListType), POINTER:: list
+    INTEGER(INTG),INTENT(OUT) :: key, value
+    INTEGER(INTG), INTENT(OUT)        :: err !<The error code.
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+
+    !Local Variables
+    TYPE(HashLinkedListItemType),POINTER :: next
+
+    ENTERS("HashLinkedList_RemoveFirst", err, error, *999)
+    
+    IF(ASSOCIATED(list%root)) THEN
+      key   = list%root%key
+      value = list%root%value
+      next => list%root%next
+      DEALLOCATE(list%root)
+      list%root => next
+      IF(ASSOCIATED(list%root)) THEN
+        IF(.NOT.ASSOCIATED(list%root%next)) list%last => list%root  ! only one left
+      ELSE
+        list%last => NULL()
+      END IF
+    ELSE
+      IF(diagnostics1) THEN
+      CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE, &
+        & ">>> warning: Linked list is empty and first item cannot be removed.",err,error,*999) 
+      ! write(*,*) ">>> warning: Linked list is empty and first item cannot be removed."
+      END IF
+    END IF
+
+
+    EXITS("HashLinkedList_RemoveFirst")
+    RETURN
+999 ERRORSEXITS("HashLinkedList_RemoveFirst",err,error)
+    RETURN
+
+  END SUBROUTINE HashLinkedList_RemoveFirst
+
+!==================================================
+
+  !> removes the last item from list and returns its value in data
+  SUBROUTINE HashLinkedList_RemoveLast(list,key,value,err,error,*)
+    
+    ! argument variables
+    TYPE(HashLinkedListType), POINTER :: list
+    INTEGER(INTG),INTENT(OUT) :: key,value
+    INTEGER(INTG), INTENT(OUT)        :: err !<The error code.
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+
+    ! local variables
+    TYPE(HashLinkedListItemType),POINTER :: current => NULL()
+
+    ENTERS("HashLinkedList_RemoveLast", err, error, *999)
+
+    IF(.NOT.ASSOCIATED(list%root)) THEN
+      IF(diagnostics1) THEN
+        CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE, &
+          & ">>> warning: Linked list is empty and last item cannot be removed.",err,error,*999) 
+      END IF
+      !write(*,*) ">>> warning: Linked list is empty and last item cannot be removed."
+      GOTO 987
+    END IF
+    current => list%root
+
+    DO
+      IF(ASSOCIATED(current%next)) THEN
+        IF(ASSOCIATED(current%next%next)) THEN
+          current => current%next
+        ELSE
+          ! next one is the last one
+          key = current%next%key
+          value = current%next%value
+          DEALLOCATE(current%next)
+          current%next => NULL()
+          list%last => current
+          EXIT
+        END IF
+      ELSE
+        ! there must be only one item in the list(?)!
+        key   = current%key
+        value = current%value
+        DEALLOCATE(list%root)
+        list%root => NULL()
+        list%last => NULL()
+        EXIT
+      END IF
+    END DO
+
+987 EXITS("HashLinkedList_RemoveLast")
+    RETURN
+999 ERRORSEXITS("HashLinkedList_RemoveLast",err,error)
+    RETURN
+
+  END SUBROUTINE HashLinkedList_RemoveLast
+
+!==================================================
+
+  !> will delete and deallocate all items
+  SUBROUTINE HashLinkedList_Destroy(list,err,error, *)
+
+    ! argument variables
+    TYPE(HashLinkedListType), POINTER :: list
+    INTEGER(INTG), INTENT(OUT)        :: err !<The error code.
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    ! local variables
+    TYPE(HashLinkedListItemType), POINTER :: current =>NULL()
+    TYPE(HashLinkedListItemType), POINTER :: next=>NULL()
+
+    ENTERS("HashLinkedList_Destroy", err, error, *999)
+
+    IF(.NOT.ASSOCIATED(list%root)) GOTO 987
+
+    current => list%root
+    DO
+      IF(ASSOCIATED(current%next)) THEN
+        next => current%next
+        DEALLOCATE(current)
+        current => next
+      ELSE
+        DEALLOCATE(current)
+        EXIT
+      END IF
+    END DO
+    
+    list%root => NULL()
+    list%last => NULL()
+    
+    ! deallocate memory
+    DEALLOCATE(list)
+
+987 EXITS("HashLinkedList_Destroy")
+    RETURN
+999 ERRORSEXITS("HashLinkedList_Destroy",err,error)
+    RETURN
+
+  END SUBROUTINE HashLinkedList_Destroy
+
+!==================================================
+
+  !> returns true if the list is empty
+  SUBROUTINE HashLinkedList_IsEmpty(list, isEmpty, err, error, *)
+
+    ! Argument Variables
+    TYPE(HashLinkedListType), POINTER :: list
+    LOGICAL, INTENT(OUT)              :: isEmpty
+    INTEGER(INTG), INTENT(OUT)        :: err !<The error code.
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+
+    ENTERS("HashLinkedList_IsEmpty", err, error, *999)
+
+    isEmpty = .TRUE.
+
+    IF(ASSOCIATED(list)) THEN
+      IF(ASSOCIATED(list%root)) isEmpty = .FALSE.
+    END IF
+
+    EXITS("HashLinkedList_IsEmpty")
+    RETURN
+999 ERRORSEXITS("HashLinkedList_IsEmpty",err,error)
+    RETURN
+
+  END SUBROUTINE HashLinkedList_IsEmpty
+
+!==================================================
+
+  !> copies out the data to an allocatable array of size 2xn
+  SUBROUTINE HashLinkedList_ListToArray(list,array,err,error,*)
+
+    ! argument variables
+    TYPE(HashLinkedListType),POINTER:: list
+
+    INTEGER(INTG),ALLOCATABLE,INTENT(OUT) :: array(:,:)
+    INTEGER(INTG), INTENT(OUT)        :: err !<The error code.
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+
+    ! local variables
+    INTEGER(INTG) :: i,n
+    LOGICAL       :: isEmpty
+    TYPE(HashLinkedListItemType),POINTER :: current => NULL()
+
+    ENTERS("HashLinkedList_ListToArray", err, error, *999)
+
+    ! return zero-size array if list is empty
+    CALL HashLinkedList_IsEmpty(list, isEmpty, err, error, *999)
+    IF(isEmpty) THEN
+      ALLOCATE(array(0,0), STAT=err)
+      IF(err/=0) CALL FlagError("Could not allocate list array",err,error,*999)
+      GOTO 987
+    END IF
+
+    ! first traversing to find size
+    current => list%root
+    n=1
+    DO
+      IF(ASSOCIATED(current%next)) THEN
+        n=n+1
+        current => current%next
+      ELSE
+        EXIT
+      END IF
+    END DO
+
+    ! copy to array
+    IF(ALLOCATEd(array)) DEALLOCATE(array)
+    ALLOCATE(array(2,n),stat=err)
+    IF(err/=0) CALL FlagError("Could not allocate list array",err,error,*999)
+    current => list%root
+    DO i=1,n
+      array(1,i)=current%key
+      array(2,i)=current%value
+
+      current => current%next
+    END DO
+
+987 EXITS("HashLinkedList_ListToArray")
+    RETURN
+999 ERRORSEXITS("HashLinkedList_ListToArray",err,error)
+    RETURN
+
+  END SUBROUTINE HashLinkedList_ListToArray
+
+!==================================================
+
+  !>Starts the creation of a list array and returns a pointer to the created list
+  SUBROUTINE HashListArray_CreateStart(listArray,err,error,*)
+
+    !Argument Variables 
+    TYPE (HashListArrayType), POINTER :: listArray !<A pointer to the list array to create. Must not be associated on entry.
+    INTEGER(INTG), INTENT(OUT) :: err !<The error code.
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    !Local Variables
+    INTEGER(INTG)        :: dummyErr !<The error code.
+    TYPE(VARYING_STRING) :: dummyError !<The error string.
+
+    ENTERS("HashListArray_CreateStart",err,error,*999)
+
+    IF(ASSOCIATED(listArray)) THEN
+      CALL FlagError("List Array is already associated",err,error,*998)
+    ELSE
+      ALLOCATE(listArray,STAT=err)
+      IF(err/=0) CALL FlagError("Could not allocate list array",err,error,*999)
+    ! CALL HashListArrayType_Initialise(tbl,err,error,*999)
+      !Set Defaults
+      listArray%vecLen=10 ! default number of elements in array
+      listArray%isFinished = .FALSE.
+    END IF
+
+    EXITS("HashListArray_CreateStart")
+    RETURN
+999 CALL HashListArray_Finalise(listArray,dummyErr,dummyError,*998)
+998 ERRORSEXITS("HashListArray_CreateStart",err,error)
+    RETURN 1
+
+  END SUBROUTINE HashListArray_CreateStart
+
+!==================================================
+
+  !>Finishes the creation of a List Array created with HashListArray_CreateStart 
+  SUBROUTINE HashListArray_CreateFinish(listArray,err,error,*)
+
+    !Argument Variables 
+    TYPE (HashListArrayType), POINTER :: listArray !<A pointer to the list array to create. Must not be associated on entry.
+    INTEGER(INTG), INTENT(OUT) :: err !<The error code.
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    !Local Variables
+    INTEGER(INTG) :: i
+
+    ENTERS("HashListArray_CreateFinish",err,error,*999)
+
+    IF(ASSOCIATED(listArray)) THEN
+       IF(listArray%isFinished) THEN
+         CALL FlagError("List Array is already finished",err,error,*999)
+       ELSE  
+         ! Allocate lists
+         ALLOCATE(listArray%vec(0:listArray%vecLen-1), STAT=err)
+         IF(err/=0) CALL FlagError("Could not allocate list array",err,error,*999)
+         ! Nullify each list pointer
+         DO i = 0, listArray%vecLen-1
+            NULLIFY(listArray%vec(i)%ptr)
+         END DO
+
+         listArray%isFinished = .TRUE.
+       END IF
+    ELSE
+       CALL FlagError("List Array is not associated",err,error,*999)
+    END IF
+
+    EXITS("HashListArray_CreateFinish")
+    RETURN
+999 ERRORSEXITS("HashListArray_CreateFinish",err,error)
+    RETURN 1
+
+  END SUBROUTINE HashListArray_CreateFinish
+
+!==================================================
+
+  !>Finalises a list array and deallocates all memory.
+  SUBROUTINE HashListArray_Finalise(listArray,err,error,*)
+
+    !Argument Variables
+    TYPE (HashListArrayType), POINTER :: listArray !<The list array to finalise
+    INTEGER(INTG), INTENT(OUT)    :: err !<The error code
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string
+
+    ENTERS("HashListArray_Finalise",err,error,*999)
+ 
+   IF(ASSOCIATED(listArray)) THEN
+      IF(ALLOCATED(listArray%vec)) DEALLOCATE(listArray%vec)
+      DEALLOCATE(listArray)
+   END IF
+
+    EXITS("HashListArray_Finalise")
+    RETURN
+999 ERRORSEXITS("HashListArray_Finalise",err,error)
+    RETURN 1
+  END SUBROUTINE HashListArray_Finalise
+
+
+  !>Finishes the creation of a table created with HashTable_CreateStart
+  SUBROUTINE HashTable_CreateFinish(table,err,error,*)
+
+    !Argument Variables
+    TYPE(HashTableType), POINTER :: table !<A pointer to the table to finish
+    INTEGER(INTG), INTENT(OUT) :: err !<The error code
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    !Local Variables
+
+!    not needed    
+!    TYPE(VARYING_STRING) :: dummyError !<The error string.
+!    INTEGER(INTG)        :: dummyErr
+
+    ENTERS("HashTable_CreateFinish",err,error,*998)
+
+    IF(ASSOCIATED(table)) THEN
+      IF(table%hashTableFinished) THEN
+        CALL FlagError("Table is already finished",err,error,*998)
+      ELSE
+         ! Allocate pointers
+         ! No pointers to allocate
+         ! Finish the tree creation
+        table%hashTableFinished=.TRUE.
+      END IF
+    ELSE
+      CALL FlagError("Table is not associated",err,error,*998)
+    END IF
+
+    EXITS("HashTable_CreateFinish")
+    RETURN
+!999 CALL HashTable_Finalise(table,dummyErr,dummyError,*998)
+! not needed
+998 ERRORSEXITS("HashTable_CreateFinish",err,error)
+    RETURN 1
+  END SUBROUTINE HashTable_CreateFinish
+
+!==================================================
+
+  !>Starts the creation of a hash table and returns a pointer to the created table.
+  SUBROUTINE HashTable_CreateStart(table,err,error,*)
+
+    !Argument Variables
+    TYPE(HashTableType), POINTER :: table !<A pointer to the table to create. Must NOT be associated on entry.
+    INTEGER(INTG), INTENT(OUT) :: err !<The error code.
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    !Local Variable
+
+    TYPE(VARYING_STRING) :: dummyError !<The error string.
+    INTEGER(INTG)        :: dummyErr
+
+    ENTERS("HashTable_CreateStart",err,error,*998)
+
+    IF(ASSOCIATED(table)) THEN
+      CALL FlagError("Table is already associated",err,error,*998)
+    ELSE
+      ALLOCATE(table,STAT=err)
+      IF(err/=0) CALL FlagError("Could not allocate table",err,error,*999)
+      CALL HashTable_Initialise(table,err,error,*999)
+    END IF
+
+    EXITS("HashTable_CreateStart")
+    RETURN
+999 CALL HashTable_Finalise(table,dummyErr,dummyError,*998)
+998 ERRORSEXITS("HashTable_CreateStart",err,error)
+    RETURN 1
+  END SUBROUTINE HashTable_CreateStart
+
+!==================================================
+
+ !>Initialises a hash table
+  SUBROUTINE HashTable_Initialise(table,err,error,*)
+
+    !Argument Variables
+    TYPE(HashTableType), POINTER :: table !<A pointer to the table to initialise
+    INTEGER(INTG), INTENT(OUT) :: err !<The error code
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    !Local Variables
+    
+    ENTERS("HashTable_Initialise",err,error,*999)
+
+    IF(ASSOCIATED(table)) THEN
+      table%hashTableFinished=.FALSE.
+      table%n=0
+      table%p=0
+      table%listSVal=>NULL()
+      IF(ALLOCATED(table%vecTKey)) DEALLOCATE(table%vecTKey)
+      IF(ALLOCATED(table%vecTVal)) DEALLOCATE(table%vecTVal)
+      IF(ALLOCATED(table%vecSKey)) DEALLOCATE(table%vecSKey)
+      IF(ALLOCATED(table%arrayOfListSVal)) DEALLOCATE(table%arrayOfListSVal)
+    ELSE
+      CALL FlagError("Table is not associated",err,error,*999)
+    END IF
+
+    EXITS("HashTable_Initialise")
+    RETURN
+999 ERRORSEXITS("HashTable_Initialise",err,error)
+    RETURN 1
+  END SUBROUTINE HashTable_Initialise
+
+!==================================================
+
+! Queries for a value q and returns index
+  SUBROUTINE HashTable_GetKey(table, q, indexFound, isFound, err, error, *) 
+
+    !Argument variables
+    TYPE(HashTableType), POINTER :: table
+    INTEGER(INTG), INTENT(IN)    :: q          ! queried value
+    LOGICAL, INTENT(OUT)         :: isFound    ! if object has been found
+    INTEGER(INTG), INTENT(OUT)   :: indexFound ! and corresponding index
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ! Local variables    
+    INTEGER(INTG) :: j, blockInd
+
+    ENTERS("HashTable_GetKey",err,error,*999)
+
+    isFound = .FALSE.
+    indexFound = 0 !Returns meaningless value (0) if q not found: Always consider flag found
+    ! copy procedure as in paper
+    j = MyAlgModule(MyAlgModule(table%vecTKey(0)*q,table%p),table%n)
+
+    blockInd = table%vecTKey(j)
+    IF(blockInd /= 0) THEN ! if block is non empty
+       indexFound = MyAlgModule(MyAlgModule(table%vecTKey(blockInd+1)*q,table%p),table%vecTKey(blockInd)**2)+1+blockInd
+       IF(table%vecTKey(indexFound)==q) isFound = .TRUE.
+    END IF
+
+    EXITS("HashTable_GetKey")
+    RETURN
+999 ERRORSEXITS("HashTable_GetKey",err,error)
+    RETURN 1
+
+  END SUBROUTINE HashTable_GetKey
+
+!==================================================
+
+  ! Returns the INT value(s) from the array of lists at key location indexFound
+  ! The list number(s) are stored in the array listNum
+  SUBROUTINE HashTable_GetListArrayIntValue(table, indexFound, value, listNum, err,error,*) 
+
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+
+    INTEGER(INTG), INTENT(IN)  :: indexFound 
+    INTEGER(INTG), ALLOCATABLE, INTENT(OUT) :: value(:,:) ! Allow for different sizes
+
+    INTEGER(INTG), ALLOCATABLE, INTENT(OUT) :: listNum(:) ! The lists where the values are found
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ! Local variables
+    INTEGER(INTG) :: indexHash
+    INTEGER(INTG), ALLOCATABLE :: listData(:)
+    INTEGER(INTG) :: listTot, listAssociated, i, j
+
+    ENTERS("HashTable_GetListArrayIntValue",err,error,*999)
+
+    IF(ASSOCIATED(table)) THEN
+      ! First get the index in the table
+      indexHash = table%vecTVal(indexFound)
+
+      ! Check if array of lists has been created
+      IF(ALLOCATED(table%arrayOfListSVal)) THEN
+        ! First list is the original list (listSVal). ListData has the dimension of the number of arrays.
+        IF(ASSOCIATED(table%arrayOfListSVal(1)%ptr)) THEN
+          ALLOCATE(listData(size(table%arrayOfListSVal,1)), STAT=err)
+          IF(err/=0) CALL FlagError("Could not allocate list data.",err,error,*999)  
+        END IF 
+
+        listTot = 0        ! The number of lists with INT values
+        listAssociated = 0 ! The total number of ASSOCIATED lists
+
+        DO i = 1,size(table%arrayOfListSVal,1)
+          IF(ASSOCIATED(table%arrayOfListSVal(i)%ptr)) THEN
+            IF(table%arrayOfListSVal(i)%ptr%DATA_TYPE==LIST_INTG_TYPE) THEN
+              listData(i)     = table%arrayOfListSVal(i)%ptr%DATA_DIMENSION ! store the dimension of INT list
+              listTot         = listTot +1 ! increase counter of INT lists
+            ELSE  
+              listData(i) = 0 ! append dimension 0
+            END IF
+            listAssociated = listAssociated +1 ! increase counter of ASSOCIATED lists
+          ELSE
+            EXIT ! Stop at first not-associated list
+          END IF
+        END DO
+
+        IF (listTot == 0) THEN
+          CALL FlagError("Could not find any data of the desired data type!",err,error,*999)
+        END IF
+
+        ALLOCATE(listNum(listTot), STAT=err)
+        IF(err/=0) CALL FlagError("Could not allocate listnum.",err,error,*999)  
+
+        ! value is an array #list X max dimension
+        ALLOCATE(value(listTot,maxval(listData(1:listAssociated),1)), STAT=err)
+        IF(err/=0) CALL FlagError("Could not allocate value.",err,error,*999)
+        value = 0
+
+        j = 1
+        ! Loop over associated lists
+        DO i = 1,listAssociated
+          SELECT CASE(listData(i)) ! Dimension of data in the list
+          CASE(0) ! Dimension 0 for not-INT list
+             ! Do nothing
+          CASE(1:)
+            CALL List_ItemGet(table%arrayOfListSVal(i)%ptr, indexHash, value(j,1:listData(i)),err,error,*999)
+            listNum(j) = i 
+            j = j+1
+          CASE DEFAULT
+            CALL FlagError("Invalid dimension of list data!",err,error,*999)
+          END SELECT
+        END DO
+       
+      ELSE
+        CALL FlagError("No extra values have been added",err,error,*999)
+      END IF
+    ELSE
+      CALL FlagError("Table is not associated",err,error,*999)
+    END IF
+
+    EXITS("HashTable_GetListArrayIntValue")
+    RETURN
+999 IF (ALLOCATED(listData)) DEALLOCATE(listData)
+    ERRORSEXITS("HashTable_GetListArrayIntValue",err,error)
+    RETURN 1
+
+  END SUBROUTINE HashTable_GetListArrayIntValue
+
+!==================================================
+
+  ! Returns the REAL SP value(s) from the array of lists at key location indexFound
+  ! The list number(s) are stored in the array listNum
+  SUBROUTINE HashTable_GetListArrayRealSpValue(table, indexFound, value, listNum, err,error,*) 
+
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+
+    INTEGER(INTG), INTENT(IN)  :: indexFound 
+    REAL(SP), ALLOCATABLE, INTENT(OUT) :: value(:,:) ! Allow for different sizes
+
+    INTEGER(INTG), ALLOCATABLE, INTENT(OUT) :: listNum(:) ! The lists where the values are found
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ! Local variables
+    INTEGER(INTG) :: indexHash
+    INTEGER(INTG), ALLOCATABLE :: listData(:)
+    INTEGER(INTG) :: listTot, listAssociated, i, j
+
+    ENTERS("HashTable_GetListArrayRealSpValue",err,error,*999)
+
+    IF(ASSOCIATED(table)) THEN
+      ! First get the index in the table
+      indexHash = table%vecTVal(indexFound)
+
+      ! Check if array of lists has been created
+      IF(ALLOCATED(table%arrayOfListSVal)) THEN
+        ! First list is the original list (listSVal). ListData has the dimension of the number of arrays.
+        IF(ASSOCIATED(table%arrayOfListSVal(1)%ptr)) THEN
+          ALLOCATE(listData(size(table%arrayOfListSVal,1)), STAT=err)
+          IF(err/=0) CALL FlagError("Could not allocate list data.",err,error,*999)  
+        END IF 
+
+        listTot = 0        ! The number of lists with REAL SP values
+        listAssociated = 0 ! The total number of ASSOCIATED lists
+
+        DO i = 1,size(table%arrayOfListSVal,1)
+          IF(ASSOCIATED(table%arrayOfListSVal(i)%ptr)) THEN
+            IF(table%arrayOfListSVal(i)%ptr%DATA_TYPE==LIST_SP_TYPE) THEN
+              listData(i)     = table%arrayOfListSVal(i)%ptr%DATA_DIMENSION ! store the dimension of REAL SP list
+              listTot         = listTot +1 ! increase counter of REAL SP lists
+            ELSE  
+              listData(i) = 0 ! append dimension 0
+            END IF
+            listAssociated = listAssociated +1 ! increase counter of ASSOCIATED lists
+          ELSE
+            EXIT ! Stop at first not-associated list
+          END IF
+        END DO
+
+        IF (listTot == 0) THEN
+          CALL FlagError("Could not find any data of the desired data type!",err,error,*999)
+        END IF
+
+        ALLOCATE(listNum(listTot), STAT=err)
+        IF(err/=0) CALL FlagError("Could not allocate listnum.",err,error,*999)  
+
+        ! value is an array #list X max dimension
+        ALLOCATE(value(listTot,maxval(listData(1:listAssociated),1)), STAT=err)
+        IF(err/=0) CALL FlagError("Could not allocate value.",err,error,*999)
+        value = 0
+
+        j = 1
+        ! Loop over associated lists
+        DO i = 1,listAssociated
+          SELECT CASE(listData(i)) ! Dimension of data in the list
+          CASE(0) ! Dimension 0 for not-REAL SP list
+             ! Do nothing
+          CASE(1:)
+            CALL List_ItemGet(table%arrayOfListSVal(i)%ptr, indexHash, value(j,1:listData(i)),err,error,*999)
+            listNum(j) = i 
+            j = j+1
+          CASE DEFAULT
+            CALL FlagError("Invalid dimension of list data!",err,error,*999)
+          END SELECT
+        END DO
+       
+      ELSE
+        CALL FlagError("No extra values have been added",err,error,*999)
+      END IF
+    ELSE
+      CALL FlagError("Table is not associated",err,error,*999)
+    END IF
+
+    EXITS("HashTable_GetListArrayRealSpValue")
+    RETURN
+999 IF (ALLOCATED(listData)) DEALLOCATE(listData)
+    ERRORSEXITS("HashTable_GetListArrayRealSpValue",err,error)
+    RETURN 1
+
+  END SUBROUTINE HashTable_GetListArrayRealSpValue
+
+
+!==================================================
+  ! Returns the SP value at index found
+  SUBROUTINE HashTable_GetRealSpValue(table, indexFound, value, err,error,*) 
+
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+
+    INTEGER(INTG), INTENT(IN)  :: indexFound 
+    REAL(SP), ALLOCATABLE, INTENT(OUT) :: value(:) ! Allow for different sizes
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ! Local variables
+    INTEGER(INTG) :: indexHash, listDataDimension
+    REAL(SP) :: value1D
+
+    ENTERS("HashTable_GetRealSpValue",err,error,*999)
+
+    IF(ASSOCIATED(table)) THEN
+
+         ! Use list (original list SVal) and the index stored in TVal
+         IF(ASSOCIATED(table%listSVal)) THEN
+
+           ! First get the index
+           indexHash = table%vecTVal(indexFound)
+           listDataDimension = table%listSVal%DATA_DIMENSION
+
+           ALLOCATE(value(listDataDimension), STAT=err)
+           IF(err/=0) CALL FlagError("Could not allocate value",err,error,*999)  
+
+           SELECT CASE(listDataDimension)
+           CASE(1)
+             CALL List_ItemGet(table%listSVal, indexHash, value1D,err,error,*999) 
+             value(1) = value1D
+           CASE(2:)
+             CALL List_ItemGet(table%listSVal, indexHash, value,err,error,*999)
+           CASE DEFAULT
+             CALL FlagError("Invalid dimension of list data!",err,error,*999)
+           END SELECT
+         ELSE
+           CALL FlagError("Values are not there!",err,error,*999)
+         END IF
+    ELSE
+       CALL FlagError("Table is not associated",err,error,*999)
+    END IF
+    
+    EXITS("HashTable_GetRealSpValue")
+    RETURN
+999 ERRORSEXITS("HashTable_GetRealSpValue",err,error)
+    RETURN 1
+
+  END SUBROUTINE HashTable_GetRealSpValue
+
+!==================================================
+  ! Returns the INT value at index found
+  SUBROUTINE HashTable_GetIntValue(table, indexFound, value, err,error,*) 
+
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+
+    INTEGER(INTG), INTENT(IN)  :: indexFound 
+    INTEGER(INTG), ALLOCATABLE, INTENT(OUT) :: value(:) ! Allow for different sizes
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ! Local variables
+    INTEGER(INTG)  :: indexHash, listDataDimension
+    INTEGER(INTG)  :: value1D
+
+    ENTERS("HashTable_GetRealSpValue",err,error,*999)
+
+    IF(ASSOCIATED(table)) THEN
+
+         ! Use list (original list SVal) and the index stored in TVal
+         IF(ASSOCIATED(table%listSVal)) THEN
+
+           ! First get the index
+           indexHash = table%vecTVal(indexFound)
+           listDataDimension = table%listSVal%DATA_DIMENSION
+
+           ALLOCATE(value(listDataDimension), STAT=err)
+           IF(err/=0) CALL FlagError("Could not allocate value",err,error,*999)  
+
+           SELECT CASE(listDataDimension)
+           CASE(1)
+             CALL List_ItemGet(table%listSVal, indexHash, value1D,err,error,*999) 
+             value(1) = value1D
+           CASE(2:)
+             CALL List_ItemGet(table%listSVal, indexHash, value,err,error,*999)
+           CASE DEFAULT
+             CALL FlagError("Invalid dimension of list data!",err,error,*999)
+           END SELECT
+         ELSE
+           CALL FlagError("Values are not there!",err,error,*999)
+         END IF
+    ELSE
+       CALL FlagError("Table is not associated",err,error,*999)
+    END IF
+    
+    EXITS("HashTable_GetIntValue")
+    RETURN
+999 ERRORSEXITS("HashTable_GetIntValue",err,error)
+    RETURN 1
+
+  END SUBROUTINE HashTable_GetIntValue
+
+!==================================================
+
+  ! Build the list of INT values according to input
+  SUBROUTINE HashTable_ListIntValuesSet (table, vecSVal, addElements, err,error,*)
+
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+
+    INTEGER(INTG), INTENT(IN) :: vecSVal(:,:)  ! array of values
+
+    LOGICAL, INTENT (IN) :: addElements
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ! Local variables
+    ! related to Lists
+    INTEGER(INTG) :: listSizeX, listSizeY, i
+    INTEGER(INTG) :: vecS1DVal(size(vecSVal,2)) ! in case only one row of values => 1D List
+
+    ENTERS("HashTable_ListIntValuesSet",err,error,*999)
+
+    ! The size of the input: the values correspond to the rows of the input array
+    ! Dimension of data
+    listSizeY = size(vecSVal,1) 
+    ! Number of entries (must correspond to number of keys, will be checked on insert)
+    listSizeX = size(vecSVal,2)
+
+    CALL HashTable_ListInitialise (table, LIST_INTG_TYPE, listsizeX, listSizeY, addElements, err,error,*999)
+
+    ! One-row case
+    IF (listSizeY==1) vecS1DVal = vecSVal(1,:) 
+
+    ! Add items (might be adding to an already filled list)
+    ! Resize before adding NEW items (otherwise ItemAdd would add unnecessary space)
+    IF (addElements) CALL List_Resize(table%listSVal, listSizeX, err, error, *999) 
+    DO i = 1,listSizeX
+        IF (listSizeY==1) THEN
+          CALL List_ItemAdd(table%listSVal,vecS1DVal(i), err, error, *999)
+        ELSE
+          CALL List_ItemAdd(table%listSVal,vecSVal(:,i), err, error, *999)
+        END IF
+    END DO
+
+    EXITS("HashTable_ListIntValuesSet")
+    RETURN
+999 ERRORSEXITS("HashTable_ListIntValuesSet",err,error)
+    RETURN 1 
+
+END SUBROUTINE
+
+!==================================================
+
+
+  ! Build the list of SP values according to input
+  SUBROUTINE HashTable_ListRealSpValuesSet (table, vecSVal, addElements, err,error,*)
+
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+
+    REAL(SP), INTENT(IN) :: vecSVal(:,:)  ! array of values
+
+    LOGICAL, INTENT (IN) :: addElements
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ! Local variables
+    ! related to Lists
+    INTEGER(INTG) :: listSizeX, listSizeY, i
+    REAL(SP) :: vecS1DVal(size(vecSVal,2)) ! in case only one row of values => 1D List
+
+    ENTERS("HashTable_ListRealSpValuesSet",err,error,*999)
+
+    ! The size of the input: the values correspond to the rows of the input array
+    ! Dimension of data
+    listSizeY = size(vecSVal,1) 
+    ! Number of entries (must correspond to number of keys, will be checked on insert)
+    listSizeX = size(vecSVal,2)
+
+    CALL HashTable_ListInitialise (table, LIST_SP_TYPE, listsizeX, listSizeY, addElements, err,error,*999)
+
+    ! One-row case
+    IF (listSizeY==1) vecS1DVal = vecSVal(1,:) 
+
+    ! Add items (can be adding to an already filled list)
+    ! Resize before adding NEW items (otherwise ItemAdd would add unnecessary space)
+    IF (addElements) CALL List_Resize(table%listSVal, listSizeX, err, error, *999) 
+    DO i = 1,listSizeX
+        IF (listSizeY==1) THEN
+         CALL List_ItemAdd(table%listSVal,vecS1DVal(i), err, error, *999)
+        ELSE
+         CALL List_ItemAdd(table%listSVal,vecSVal(:,i), err, error, *999)
+        END IF
+    END DO
+
+    EXITS("HashTable_ListRealSpValuesSet")
+    RETURN
+999 ERRORSEXITS("HashTable_ListRealSpValuesSet",err,error)
+    RETURN 1 
+
+END SUBROUTINE
+
+!==================================================
+
+  ! Initialise the list of values independent of data type
+  SUBROUTINE HashTable_ListInitialise (table, dataType, listsizeX, listSizeY, addElements, err,error,*)
+
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+    INTEGER(INTG), INTENT(IN)    :: dataType, listsizeX, listSizeY
+    LOGICAL, INTENT(IN)          :: addElements
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+
+    ENTERS("HashTable_ListInitialise",err,error,*999)
+
+    IF(ASSOCIATED(table)) THEN
+      IF(table%hashTableFinished) THEN
+       ! Do nothing 
+      ELSE
+        CALL FlagError("The table has not been finished.",err,error,*999)
+      END IF
+    ELSE
+      CALL FlagError("Table not associated.",err,error,*999)
+    END IF
+
+    ! Check if extra values are there
+    IF (.NOT. addElements) THEN
+      ! Deallocate extra values if creating a new table (and go on)
+      IF (ALLOCATED(table%arrayOfListSVal)) DEALLOCATE(table%arrayOfListSVal)
+    ELSE
+      IF (ALLOCATED(table%arrayOfListSVal)) THEN
+        DEALLOCATE(table%arrayOfListSVal)
+        CALL FlagError("Cannot add new elements to a list with extra values! (FIX?)",err,error,*999)
+      END IF
+    END IF
+    
+
+    IF (.NOT. addElements .AND. ASSOCIATED(table%listSVal)) CALL List_Destroy(table%listSVal, err, error, *999) ! destroy the existing list 
+    ! If not already existing, create a new one
+    IF (.NOT. ASSOCIATED(table%listSVal)) THEN
+    ! Pointers nullified in hash initialise!
+      CALL List_CreateStart(table%listSVal, err, error, *999)
+!   Defaults in initialise:
+!      LIST%LIST_FINISHED=.FALSE.
+!      LIST%MUTABLE=.FALSE.
+!      LIST%NUMBER_IN_LIST=0
+!      LIST%DATA_DIMENSION=1
+!      LIST%INITIAL_SIZE=10
+!      LIST%SIZE=0
+!      LIST%DATA_TYPE=LIST_INTG_TYPE
+!      LIST%KEY_DIMENSION=1
+!      LIST%SORT_ORDER=LIST_SORT_ASCENDING_TYPE
+!      LIST%SORT_METHOD=LIST_HEAP_SORT_METHOD
+
+    ! Let us create a list of intg of dim listSizeY
+      CALL List_DataDimensionSet(table%listSVal,listSizeY, err, error, *999)
+    ! Initial size listSizeX
+      CALL List_InitialSizeSet(table%listSVal, listSizeX, err, error, *999)
+      CALL List_DataTypeSet(table%listSVal, dataType, err, error, *999) ! Ok default
+      CALL List_CreateFinish(table%listSVal, err, error, *999)
+    END IF
+
+    EXITS("HashTable_ListInitialise")
+    RETURN
+999 ERRORSEXITS("HashTable_ListInitialise",err,error)
+    RETURN 1 
+
+END SUBROUTINE
+
+!==================================================
+
+  ! Build the list of additional INT values according to input
+  SUBROUTINE HashTable_AdditionalIntValuesSet (table, vecSVal, dataType, err,error,*)
+
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+
+    INTEGER(INTG), INTENT(IN) :: vecSVal(:,:)  ! array of values
+    INTEGER(INTG) :: dataType  ! data type
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ! Local variables
+    ! related to Lists
+    INTEGER(INTG) :: listSizeX, listSizeY, listNum, j
+    INTEGER(INTG) :: vecS1DVal(size(vecSVal,2)) ! in case only one row of values => 1D List
+
+    ENTERS("HashTable_AdditionalIntValuesSet",err,error,*999)
+
+    ! The size of the input: the values correspond to the rows of the input array
+    ! Dimension of data
+    listSizeY = size(vecSVal,1) 
+    ! Number of entries (must correspond to number of keys, will be checked on insert)
+    listSizeX = size(vecSVal,2)
+
+    ! Init new list and return number in the list array
+    CALL HashTable_AdditionalValuesInit (table, dataType, listSizeX, listSizeY, listNum, err,error,*999)
+
+    IF (listSizeY==1) vecS1DVal = vecSVal(1,:) 
+
+    ! Add items 
+    DO j = 1,listSizeX
+      IF (listSizeY==1) THEN
+        CALL List_ItemAdd(table%arrayOfListSVal(listNum)%ptr,vecS1DVal(j), err, error, *999)
+      ELSE
+        CALL List_ItemAdd(table%arrayOfListSVal(listNum)%ptr,vecSVal(:,j), err, error, *999)
+      END IF
+    END DO
+
+    EXITS("HashTable_AdditionalIntValuesSet")
+    RETURN
+999 ERRORSEXITS("HashTable_AdditionalIntValuesSet",err,error)
+    RETURN 1 
+
+END SUBROUTINE
+
+!==================================================
+
+  ! Build the list of additional REAL SP values according to input
+  SUBROUTINE HashTable_AdditionalRealSpValuesSet (table, vecSVal, dataType, err,error,*)
+
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+
+    REAL(SP), INTENT(IN) :: vecSVal(:,:)  ! array of values
+    INTEGER(INTG) :: dataType  ! data type
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ! Local variables
+    ! related to Lists
+    INTEGER(INTG) :: listSizeX, listSizeY, listNum, j
+    REAL(SP) :: vecS1DVal(size(vecSVal,2)) ! in case only one row of values => 1D List
+
+    ENTERS("HashTable_AdditionalRealSpValuesSet",err,error,*999)
+
+    ! The size of the input: the values correspond to the rows of the input array
+    ! Dimension of data
+    listSizeY = size(vecSVal,1) 
+    ! Number of entries (must correspond to number of keys, will be checked on insert)
+    listSizeX = size(vecSVal,2)
+
+    ! Init new list and return number in the list array
+    CALL HashTable_AdditionalValuesInit (table, dataType, listSizeX, listSizeY, listNum, err,error,*999)
+
+    IF (listSizeY==1) vecS1DVal = vecSVal(1,:) 
+
+    ! Add items 
+    DO j = 1,listSizeX
+      IF (listSizeY==1) THEN
+        CALL List_ItemAdd(table%arrayOfListSVal(listNum)%ptr,vecS1DVal(j), err, error, *999)
+      ELSE
+        CALL List_ItemAdd(table%arrayOfListSVal(listNum)%ptr,vecSVal(:,j), err, error, *999)
+      END IF
+    END DO
+
+    EXITS("HashTable_AdditionalRealSpValuesSet")
+    RETURN
+999 ERRORSEXITS("HashTable_AdditionalRealSpValuesSet",err,error)
+    RETURN 1 
+
+END SUBROUTINE
+
+!==================================================
+
+  ! Initialise additional list of values (returns newly allocated list)
+  SUBROUTINE HashTable_AdditionalValuesInit (table, dataType, listSizeX, listSizeY, listNum, err,error,*)
+
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+
+    INTEGER(INTG), INTENT(IN)  :: dataType, listSizeX, listSizeY
+    INTEGER(INTG), INTENT(OUT) :: listNum
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ! Local variables
+    INTEGER(INTG) :: i, sizeInit
+
+    ENTERS("HashTable_AdditionalValuesInit",err,error,*999)
+   
+   ! New array of lists if none allocated yet
+    IF (.NOT. ALLOCATED(table%arrayOfListSVal)) THEN
+
+      ! Allocate a sufficient dimension of (empty) lists
+      ALLOCATE(table%arrayOfListSVal(HASH_INITIAL_ADDITIONAL_VALUES_SIZE),STAT=err)   
+      DO i=1,HASH_INITIAL_ADDITIONAL_VALUES_SIZE
+        table%arrayOfListSVal(i)%ptr=>NULL() ! Define the ptr to list as not associated
+      END DO
+    
+      CALL List_NumberOfItemsGet(table%listSVal, sizeInit, err, error, *999)
+ 
+      IF (sizeInit /= listSizeX) THEN
+        DEALLOCATE(table%arrayOfListSVal)
+        CALL FlagError("Adding values of different size than keys in table!!", err, error,*999)
+      END IF
+
+      ! First list in array is the one already provided in listSval (at table creation)
+      ! No need to initialise / finish new pointer. Just allocate.
+      ALLOCATE(table%arrayOfListSVal(1)%ptr,STAT=err)
+      IF(err/=0) CALL FlagError("Could not allocate array list 1.", err, error,*999)
+      table%arrayOfListSVal(1)%ptr = table%listSval 
+
+    END IF
+
+    ! Allocate an additional list (look for the first free spot)
+    listNum = 2
+    DO
+      IF (size(table%arrayOfListSVal)>=listNum) THEN
+        IF (.NOT. ASSOCIATED(table%arrayOfListSVal(listNum)%ptr)) THEN
+          CALL List_CreateStart(table%arrayOfListSVal(listNum)%ptr, err, error, *999)
+!   Defaults in initialise:
+!      LIST%LIST_FINISHED=.FALSE.
+!      LIST%MUTABLE=.FALSE.
+!      LIST%NUMBER_IN_LIST=0
+!      LIST%DATA_DIMENSION=1
+!      LIST%INITIAL_SIZE=10
+!      LIST%SIZE=0
+!      LIST%DATA_TYPE=LIST_INTG_TYPE
+!      LIST%KEY_DIMENSION=1
+!      LIST%SORT_ORDER=LIST_SORT_ASCENDING_TYPE
+!      LIST%SORT_METHOD=LIST_HEAP_SORT_METHOD
+
+       ! Let us create a list of intg of dim listSizeY
+          CALL List_DataDimensionSet(table%arrayOfListSVal(listNum)%ptr, listSizeY, err, error, *999)
+          CALL List_DataTypeSet     (table%arrayOfListSVal(listNum)%ptr, dataType, err, error, *999) 
+       ! Initial size listSizeX
+          CALL List_InitialSizeSet  (table%arrayOfListSVal(listNum)%ptr, listSizeX, err, error, *999)
+          CALL List_CreateFinish    (table%arrayOfListSVal(listNum)%ptr, err, error, *999)
+          EXIT
+        END IF
+        listNum = listNum+1
+      ELSE
+        CALL FlagError("Number of additional lists exceeds initial allocation!",err,error,*999)
+      END IF      
+    END DO
+
+    EXITS("HashTable_AdditionalValuesInit")
+    RETURN
+999 ERRORSEXITS("HashTable_AdditionalValuesInit",err,error)
+    RETURN 1 
+
+END SUBROUTINE
+
+!==================================================
+
+    ! Create the list of values and call the hash table create routine 
+    SUBROUTINE HashTable_IntValuesSetAndInsert(table,  vecSKey, vecSVal, addElements, err,error,*)
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+
+    INTEGER(INTG), INTENT(IN) :: vecSKey(:)    ! array of keys
+    INTEGER(INTG), INTENT(IN) :: vecSVal(:,:)  ! array of values
+
+    LOGICAL, INTENT(IN) :: addElements ! adding to existing or new table
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ENTERS("HashTable_IntMultValuesSetAndInsert",err,error,*999)
+
+    ! Create the list of values (if empty)
+    CALL HashTable_ListValuesSet(table, vecSVal, addElements, err,error,*999)
+    ! Compute the table based on the key array
+    CALL HashTable_Insert(table, vecSKey, addElements, err,error,*999)
+
+    EXITS("HashTable_IntMultValuesSetAndInsert")
+    RETURN
+999 ERRORSEXITS("HashTable_IntMultValuesSetAndInsert",err,error)
+    RETURN 1 
+
+    END SUBROUTINE
+
+!==================================================
+
+    ! Create the list of values and call the hash table create routine 
+    SUBROUTINE HashTable_RealSpValuesSetAndInsert(table,  vecSKey, vecSVal, addElements, err,error,*)
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+
+    INTEGER(INTG), INTENT(IN) :: vecSKey (:)  ! array of keys
+    REAL(SP), INTENT(IN) :: vecSVal(:,:)      ! array of values
+
+    LOGICAL, INTENT(IN) :: addElements ! adding to existing or new table
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ENTERS("HashTable_RealSpMultValuesSetAndInsert",err,error,*999)
+
+    ! Create the list of values (if empty)
+    CALL HashTable_ListValuesSet(table, vecSVal, addElements, err,error,*999)
+    ! Compute the table based on the key array
+    CALL HashTable_Insert(table, vecSKey, addElements, err,error,*999)
+
+    EXITS("HashTable_RealSpMultValuesSetAndInsert")
+    RETURN
+999 ERRORSEXITS("HashTable_RealSpMultValuesSetAndInsert",err,error)
+    RETURN 1 
+
+    END SUBROUTINE
+
+
+!==================================================
+
+
+  ! Build the table according to Fredman 1984: This is the most important subroutine
+   SUBROUTINE HashTable_Insert(table, vecSKey, addElements, err,error,*)
+    ! Argument variables
+    TYPE(HashTableType), POINTER :: table    
+
+    INTEGER(INTG), INTENT(IN) :: vecSKey(:)    ! array of keys
+
+    LOGICAL, INTENT(IN) :: addElements ! adding to existing or new table
+
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+    INTEGER(INTG), INTENT(OUT) :: err
+
+    ! Local variables
+    ! related to subset W
+    TYPE (HashListArrayType), POINTER  :: subsetW=>NULL() ! W as a pointer to an array of lists
+    TYPE(HashLinkedListItemType) :: subsetWListItem
+    INTEGER(INTG) :: cardSubsetW, storeCardSubsetW, checkSumW, checkSquaredSumW, countFilledSubsetW, indexSubsetW
+    ! related to array T
+    INTEGER(INTG) :: indexTLoop, sizeOfT, prevIndexT
+    INTEGER(INTG), ALLOCATABLE :: indexTArray (:), tmpKey(:)
+    ! related to Lists
+    INTEGER(INTG) :: listSize
+
+    TYPE(VARYING_STRING) :: dummyError !<The error string.
+    INTEGER(INTG)        :: dummyErr
+    ! generic
+    INTEGER(INTG) :: i,k,j
+    LOGICAL :: isFound = .FALSE.
+
+    ENTERS("HashTable_Insert",err,error,*999)
+
+    IF(addElements) THEN ! add to existing table
+     
+      IF(.NOT. ASSOCIATED(table))        CALL FlagError("Table should be already created.",err,error,*999)
+      IF(.NOT. table%hashTableFinished)  CALL FlagError("Table should be already initialised.",err,error,*999) 
+      IF(ALLOCATED(table%vecTKey)) DEALLOCATE(table%vecTKey)
+      IF(ALLOCATED(table%vecTVal)) DEALLOCATE(table%vecTVal)
+
+      IF(.NOT. ALLOCATED(table%vecSKey)) CALL FlagError("Cannot put new keys in an empty table.",err,error,*999)
+      IF(.NOT. ASSOCIATED(table%listSVal)) CALL FlagError("Cannot put new values in an empty table.",err,error,*999)
+
+      table%n = size(table%vecSKey,1)  ! old size of table
+
+      ALLOCATE(tmpKey(table%n+size(vecSKey,1)), STAT=err)
+      IF(err/=0) CALL FlagError("Could not allocate tmp in table.",err,error,*999)
+      tmpKey(1:table%n)  = table%vecSKey
+      tmpKey(table%n+1:) = vecSKey(1:)
+      CALL MOVE_ALLOC (tmpKey, table%vecSKey)
+
+      table%n = size(table%vecSKey,1) ! new size of table
+
+    ELSE ! new table
+
+      IF(ASSOCIATED(table)) THEN
+        IF(table%hashTableFinished) THEN
+      !   DEALLOCATE already existing memory in table
+          IF(ALLOCATED(table%vecTKey)) DEALLOCATE(table%vecTKey)
+          IF(ALLOCATED(table%vecTVal)) DEALLOCATE(table%vecTVal)
+          IF(ALLOCATED(table%vecSKey)) DEALLOCATE(table%vecSKey)
+        ELSE
+          CALL FlagError("The table has not been finished.",err,error,*999)
+        END IF
+      ELSE
+        CALL FlagError("Table not associated.",err,error,*999)
+      END IF
+
+      table%n = size(vecSKey,1) 
+
+      ALLOCATE(table%vecSKey(table%n), STAT=err)
+      IF(err/=0) CALL FlagError("Could not allocate vecSKey in table",err,error,*999)
+      table%vecSKey(1:table%n) = vecSKey(1:table%n)
+
+    END IF
+
+
+    ! Check that size of list corresponds to key size
+     IF (ASSOCIATED(table%listSVal)) THEN
+       CALL List_NumberOfItemsGet(table%listSVal, listSize,err, error, *999)
+       IF(table%n /= listSize)  CALL FlagError("Keys and list of values have different size!",err,error,*999)
+     END IF       
+  
+
+    ! Compute p as the first prime number larger than biggest key
+    table%p = maxval(table%vecSKey,1)+1
+    SELECT CASE(table%p)
+    CASE (1)
+    CASE (2)
+    CASE DEFAULT
+     i=2
+     DO !i = 2, table%p-1
+        IF(MOD(table%p,i) == 0)  THEN
+          table%p = table%p+1
+          i = 2
+          CYCLE
+        END IF
+        i = i+1
+        IF(i == table%p) EXIT
+     END DO
+    END SELECT 
+
+    ! Initialise w as an array of size n
+    CALL HashListArray_CreateStart(subsetW, err,error,*999)
+    subsetW%vecLen = table%n
+    CALL HashListArray_CreateFinish(subsetW, err,error,*999)
+
+    ! Compute the W_j
+    k=2 ! use 2 as start as in paper (but 1 also works for paper case!!)
+
+    DO
+
+     DO i=1,table%n
+
+      indexSubsetW = MyAlgModule(MyAlgModule((k*table%vecSKey(i)),table%p),table%n)
+      !PRINT *, indexSubsetW
+      ! CALL HashLinkedList_AddData(subsetW%vec(indexSubsetW-1)%ptr,table%vecSKey(i),table%vecSVal(i), isFound, err, error, *999) ! numbered from 0!
+      ! Change! Store INDEX i instead of value!
+      CALL HashLinkedList_AddData(subsetW%vec(indexSubsetW-1)%ptr,table%vecSKey(i), i, isFound, err, error, *999) ! numbered from 0!
+      IF(isFound) CALL FlagError("Duplicate key, error!",err,error,*999) 
+         !PRINT *, "Success with ", indexSubsetW
+     END DO
+
+     checkSquaredSumW =0
+     checkSumW   = 0
+     countFilledSubsetW  = 0
+     DO i=1,table%n
+      CALL HashLinkedList_Size(subsetW%vec(i-1)%ptr, cardSubsetW, err, error,*999) 
+      IF(cardSubsetW /= 0) THEN
+      checkSquaredSumW = checkSquaredSumW + cardSubsetW**2
+      checkSumW    = checkSumW    + cardSubsetW
+      countFilledSubsetW   = countFilledSubsetW +1
+      END IF
+     END DO
+
+    IF(checkSquaredSumW<3*table%n) THEN ! k is correct
+     EXIT
+    ELSE
+     k = k+1 ! try another k
+     IF(k>table%p) CALL FlagError("Suitable k could not be found!",err,error,*999) 
+    END IF
+
+   END DO
+
+   ! Allocate T: T(0)=k, T(1),...,T(n) access indices, T(*) |W_j| T(**) k' T(***) = W_j
+   sizeOfT = 1+table%n+2*countFilledSubsetW+checkSquaredSumW
+   ALLOCATE (table%vecTKey(0:(sizeOfT-1)),STAT=err)
+   IF(err/=0) CALL FlagError("Could not allocate vecTKey in table.",err,error,*999)
+   ALLOCATE (table%vecTVal(0:(sizeOfT-1)),STAT=err)
+   IF(err/=0) CALL FlagError("Could not allocate vecTVal in table.",err,error,*999)
+   table%vecTKey = 0
+   table%vecTVal = 0
+   table%vecTKey(0) = k
+   IF(sizeOfT<=table%n+1) CALL FlagError("No keys to be assigned!",err,error,*999) 
+
+   prevIndexT = table%n
+
+   DO i=1,table%n
+
+     ! Cardinality of block W_i 
+     CALL HashLinkedList_Size(subsetW%vec(i-1)%ptr, storeCardSubsetW, err,error,*999)
+
+     IF(storeCardSubsetW == 0) THEN
+       table%vecTKey(i) = 0 ! Block is empty
+     ELSE
+
+       IF(prevIndexT == table%n) THEN
+        table%vecTKey(i) = table%n+1 ! First non-zero entry is just index that follows n
+       ELSE 
+        ! Compute card of 
+        ! CALL HashLinkedList_Size(subsetW%vec(prevIndexT-1), cardSubsetW, err)    
+        ! replaced with cardSubsetW = storeCardSubsetW below!!!!   
+        table%vecTKey(i) = table%vecTKey(prevIndexT)+2+cardSubsetW**2    
+       END IF
+
+       table%vecTKey(table%vecTKey(i))   = storeCardSubsetW ! First  element of each block: card of W
+       table%vecTKey(table%vecTKey(i)+1) = 1            ! Second element of each block: k' (set to 1)
+
+       ALLOCATE (indexTArray(storeCardSubsetW), STAT=err)  ! indexTArray contains the indices where keys will be stored in T for each block    
+       IF(err/=0) CALL FlagError("Could not allocate indexTArray in table.",err,error,*999)
+       subsetWListItem = subsetW%vec(i-1)%ptr%root            ! take the first item in the block
+       j = 1 ! init j for do loop
+
+       DO 
+         ! function given by corollary 2
+         ! gives LOCAL index
+         ! +2 translation
+         indexTArray(j) =  MyAlgModule(MyAlgModule(table%vecTKey(table%vecTKey(i)+1)*subsetWListItem%key,table%p) &
+                             & ,storeCardSubsetW**2)+2
+         IF((j==1).OR. any(indexTArray(1:(j-1))/=indexTArray(j))) THEN ! insert items until conflict
+
+            table%vecTKey(indexTArray(j)+table%vecTKey(i)-1) = subsetWListItem%key 
+            table%vecTVal(indexTArray(j)+table%vecTKey(i)-1) = subsetWListItem%value
+
+            IF(ASSOCIATED(subsetWListItem%next)) subsetWListItem = subsetWListItem%next
+            
+            IF(j==storeCardSubsetW) EXIT ! everything filled correctly
+            j = j+1  ! next key
+            
+         ELSE ! conflict: increase k' and start over
+             table%vecTKey(table%vecTKey(i)+1) =table%vecTKey(table%vecTKey(i)+1)+1 
+
+             DO indexTLoop =1,size(indexTArray,1) 
+                table%vecTKey(indexTArray(indexTLoop) + table%vecTKey(i)-1) = 0   ! reset to zero what already filled (keys)
+                table%vecTVal(indexTArray(indexTLoop) + table%vecTKey(i)-1) = 0   ! reset to zero what already filled (values)
+             END DO 
+
+             indexTArray = 0                       ! reset to zero indices
+             j = 1                              ! reset j
+             subsetWListItem = subsetW%vec(i-1)%ptr%root  ! go back to list start
+         END IF
+ 
+       END DO
+       
+       !PRINT *, "IND IN T",  indexTArray
+       DEALLOCATE(indexTArray)
+
+       prevIndexT = i
+       cardSubsetW = storeCardSubsetW
+     END IF
+
+   END DO
+
+   ! Destroy the list array of Ws
+   CALL HashListArray_Finalise(subsetW, dummyErr,dummyError,*998)
+
+   ! DEALLOCATE local memory
+   IF(ALLOCATED(indexTArray)) DEALLOCATE(indexTArray)
+
+   ! Diagnostics instead of print
+   IF(diagnostics1) THEN
+
+     CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "The number of keys n:    ", table%n, err, error, *999)
+
+     CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "The size of the array T: ", size(table%vecTKey,1), err, error, *999)
+
+     CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE, "The value of k = T(0):   ", table%vecTkey(0), err, error, *999)
+
+     CALL WriteStringVector(DIAGNOSTIC_OUTPUT_TYPE, 2,1, table%n+1, 4,4, table%vecTKey, & 
+       & '("T(1) to T(n) (indices to T blocks):", 4(x,I8))','(35x,4(x,I8))', err, error, *999)
+
+     CALL WriteStringVector(DIAGNOSTIC_OUTPUT_TYPE, table%n+2,1, size(table%vecTKey,1), 4,4, table%vecTKey, & 
+       & '("T(n+1) to T(end):                  ", 4(x,I8))','(35x,4(x,I8))', err, error, *999)
+
+!       print *, "T vector n    ", table%vecTKey 
+  !>Writes the given INTEGER VECTOR to the given output stream specified by ID. The FIRST_FORMAT is the format initially used, followed by the REPEAT_FORMAT which is repeated as many times as necessary. NUMBER_FIRST is the number of data items in the FIRST_FORMAT and NUMBER_REPEAT is the number of data items in the REPEAT_FORMAT. FIRST_IDX and LAST_IDX are the extents of the data and DELTA is the NUMBER of indices to skip for each index.
+!          CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,DOMAIN_LINE%BASIS%NUMBER_OF_NODES,4,4,DOMAIN_LINE%NODES_IN_LINE, &
+!           & '("        Nodes in line        :",4(X,I8))','(30X,4(X,I8))',err,error,*999)
+
+!  SUBROUTINE WRITE_STRING_VECTOR_DP(ID,FIRST_IDX,DELTA,LAST_IDX,                               NUMBER_FIRST,NUMBER_REPEAT,VECTOR,FIRST_FORMAT,REPEAT_FORMAT, &
+!    & err,error,*)
+
+       !print *, "n    ", table%n
+       !print *, "T vector n    ", table%vecTKey (1:table%n)
+       !print *, "T vector rest ", table%vecTKey (table%n+1:)
+       !print *, "T value vector rest ", table%vecTVal (table%n+1:)
+       !print *, "of size ", size(table%vecTKey)
+  END IF
+
+    EXITS("HashTable_Insert")
+    RETURN
+999 CALL HashTable_Finalise(Table,dummyErr,dummyError,*998)
+998 ERRORSEXITS("HashTable_Insert",err,error)
+    RETURN 1
+
+  END SUBROUTINE HashTable_Insert
+
+!==================================================
+
+  !>Finalises a hash table and deallocates all memory.
+  SUBROUTINE HashTable_Finalise(table,err,error,*)    
+
+    !Argument Variables
+    TYPE (HashTableType), POINTER :: table !<The table to finalise
+    INTEGER(INTG), INTENT(OUT)    :: err !<The error code
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string
+
+    ENTERS("HashTable_Finalise",err,error,*999)
+ 
+   IF(ASSOCIATED(table)) THEN
+      IF(ALLOCATED(table%vecTKey)) DEALLOCATE(table%vecTKey)
+      IF(ALLOCATED(table%vecTVal)) DEALLOCATE(table%vecTVal)
+      IF(ALLOCATED(table%vecSKey)) DEALLOCATE(table%vecSKey)
+      DEALLOCATE(table)
+   END IF
+
+    EXITS("HashTable_Finalise")
+    RETURN
+999 ERRORSEXITS("HashTable_Finalise",err,error)
+    RETURN 1
+
+  END SUBROUTINE HashTable_Finalise
+
+!==================================================
+
+! Computes modulo as in Fredman84 (returns divisor instead of 0)
+FUNCTION MyAlgModule(a,b) RESULT(c)
+
+    INTEGER(INTG), INTENT(IN)  :: a,b
+    INTEGER(INTG) :: c
+
+    c = MOD(a,b)
+    IF(c==0) c = b
+
+END FUNCTION MyAlgModule
+
+
+END MODULE HashRoutines
diff --git a/src/linkedlist_routines.F90 b/src/linkedlist_routines.F90
index 4f531f76..c8c6d2bc 100644
--- a/src/linkedlist_routines.F90
+++ b/src/linkedlist_routines.F90
@@ -34,6 +34,9 @@ MODULE LinkedList_routines
   public :: LinkedList_Add,LinkedList_Destroy,LinkedList_Remove_First,LinkedList_Remove_Last
   public :: LinkedList_is_Empty,LinkedList_to_Array
 
+
+
+
 contains
 
 ! -------------------------------------------------------------------
@@ -72,7 +75,7 @@ End SUBROUTINE LinkedList_Add_Data
 
 ! -------------------------------------------------------------------
 
-  !> adds all data from one list to another
+  !> adds all data from one list (addlist) to another (list)
   Subroutine LinkedList_Add_List(list,addlist,ERR,ERROR,*)
     type(LinkedList),intent(inout) :: list
     type(LinkedList),intent(in) :: addlist
@@ -157,7 +160,7 @@ Subroutine LinkedList_Remove_Last(list,data,ERR,ERROR,*)
           exit
         endif
       else
-        ! there must be only one item in the list?
+        ! there must be only one item in the list(?)!
         data = current%data
         deallocate(list%root)
         list%root => NULL()
@@ -252,4 +255,7 @@ Subroutine LinkedList_to_Array(list,array,ERR,ERROR,*)
 
   End Subroutine LinkedList_to_Array
 
+  
+
+
 End Module LinkedList_routines
diff --git a/src/lists.F90 b/src/lists.F90
index b8e1f10d..5a9131e4 100644
--- a/src/lists.F90
+++ b/src/lists.F90
@@ -112,7 +112,7 @@ MODULE LISTS
     MODULE PROCEDURE LIST_DATA_TYPE_SET
   END INTERFACE List_DataTypeSet
   
-  !>Detaches the list values from a list and returns them as a pointer to a array of base type before destroying the list \see Lists.
+  !>Detaches the list values from a list and returns them as a pointer to an array of base type before destroying the list \see Lists.
   INTERFACE LIST_DETACH_AND_DESTROY
     MODULE PROCEDURE LIST_DETACH_AND_DESTROY_INTG1
     MODULE PROCEDURE LIST_DETACH_AND_DESTROY_INTG2
@@ -122,7 +122,7 @@ MODULE LISTS
     MODULE PROCEDURE LIST_DETACH_AND_DESTROY_DP2
   END INTERFACE LIST_DETACH_AND_DESTROY
 
-  !>Detaches the list values from a list and returns them as a pointer to a array of base type before destroying the list \see Lists.
+  !>Detaches the list values from a list and returns them as a pointer to an array of base type before destroying the list \see Lists.
   INTERFACE List_DetachAndDestroy
     MODULE PROCEDURE LIST_DETACH_AND_DESTROY_INTG1
     MODULE PROCEDURE LIST_DETACH_AND_DESTROY_INTG2
@@ -412,6 +412,8 @@ MODULE LISTS
   PUBLIC LIST_SUBSET_OF
   
   PUBLIC List_SubsetOf
+ 
+  PUBLIC List_Resize
 
 CONTAINS
 
@@ -548,7 +550,7 @@ END SUBROUTINE LIST_DATA_DIMENSION_SET
   !================================================================================================================================
   !
 
-  !>Sets/changes the data dimension for a list.
+  !>Sets/changes the mutability for a list.
   SUBROUTINE LIST_MUTABLE_SET(LIST,MUTABLE,ERR,ERROR,*)
 
     !Argument Variables
@@ -898,6 +900,110 @@ SUBROUTINE LIST_INITIAL_SIZE_SET(LIST,INITIAL_SIZE,ERR,ERROR,*)
     RETURN 1
   END SUBROUTINE LIST_INITIAL_SIZE_SET
 
+  !
+  !================================================================================================================================
+  !
+
+  !>Resizes the list with a given size
+  ! Should be called BEFORE item_add to avoid doubling size 
+  SUBROUTINE List_Resize(list,newAddSize,err,error,*)
+   !Argument Variables
+    TYPE(LIST_TYPE), POINTER, INTENT(INOUT) :: LIST !<A pointer to the list
+    INTEGER(INTG), INTENT(IN) :: newAddSize !<The new size
+    INTEGER(INTG), INTENT(OUT) :: ERR !<The error code
+    TYPE(VARYING_STRING), INTENT(OUT) :: ERROR !<The error string
+    !Local Variables
+    INTEGER(INTG) :: newSize
+    INTEGER(INTG), ALLOCATABLE :: NEW_LIST_INT1(:)
+    INTEGER(INTG), ALLOCATABLE :: NEW_LIST_INT2(:,:)
+    REAL(SP), ALLOCATABLE :: NEW_LIST_SP1(:)
+    REAL(SP), ALLOCATABLE :: NEW_LIST_SP2(:,:)
+    REAL(DP), ALLOCATABLE :: NEW_LIST_DP1(:)
+    REAL(DP), ALLOCATABLE :: NEW_LIST_DP2(:,:)
+
+    TYPE(VARYING_STRING) :: LOCAL_ERROR
+
+    ENTERS("List_Resize",ERR,ERROR,*999)
+
+    IF(ASSOCIATED(LIST)) THEN
+      IF(LIST%LIST_FINISHED) THEN
+        IF(LIST%NUMBER_IN_LIST==LIST%SIZE) THEN
+          newSize = LIST%NUMBER_IN_LIST + newAddSize   
+
+          SELECT CASE (LIST%DATA_TYPE)
+          CASE (LIST_INTG_TYPE)
+            SELECT CASE (LIST%DATA_DIMENSION)
+            CASE(1)
+              ALLOCATE(NEW_LIST_INT1(newSize),STAT=ERR)
+              IF(ERR/=0) CALL FlagError("Could not allocate new list.",ERR,ERROR,*999)
+              NEW_LIST_INT1(1:LIST%NUMBER_IN_LIST)=LIST%LIST_INTG(1:LIST%NUMBER_IN_LIST)
+              CALL MOVE_ALLOC(NEW_LIST_INT1,LIST%LIST_INTG)
+            CASE(2:) 
+              ALLOCATE(NEW_LIST_INT2(LIST%DATA_DIMENSION,newSize),STAT=ERR)
+              IF(ERR/=0) CALL FlagError("Could not allocate new list.",ERR,ERROR,*999)
+              NEW_LIST_INT2(:,1:LIST%NUMBER_IN_LIST)=LIST%LIST_INTG2(:,1:LIST%NUMBER_IN_LIST)
+              CALL MOVE_ALLOC(NEW_LIST_INT2,LIST%LIST_INTG2)
+            CASE DEFAULT
+              CALL FlagError("Invalid data dimension.",ERR,ERROR,*999)
+            END SELECT
+          CASE (LIST_SP_TYPE)  
+            SELECT CASE (LIST%DATA_DIMENSION)
+            CASE(1)
+              ALLOCATE(NEW_LIST_SP1(newSize),STAT=ERR)
+              IF(ERR/=0) CALL FlagError("Could not allocate new list.",ERR,ERROR,*999)
+              NEW_LIST_SP1(1:LIST%NUMBER_IN_LIST)=LIST%LIST_SP(1:LIST%NUMBER_IN_LIST)
+              CALL MOVE_ALLOC(NEW_LIST_SP1,LIST%LIST_SP)
+            CASE(2:) 
+              ALLOCATE(NEW_LIST_SP2(LIST%DATA_DIMENSION,newSize),STAT=ERR)
+              IF(ERR/=0) CALL FlagError("Could not allocate new list.",ERR,ERROR,*999)
+              NEW_LIST_SP2(:,1:LIST%NUMBER_IN_LIST)=LIST%LIST_SP2(:,1:LIST%NUMBER_IN_LIST)
+              CALL MOVE_ALLOC(NEW_LIST_SP2,LIST%LIST_SP2)
+            CASE DEFAULT
+              CALL FlagError("Invalid data dimension.",ERR,ERROR,*999)
+            END SELECT
+          CASE (LIST_DP_TYPE)
+            SELECT CASE (LIST%DATA_DIMENSION)
+            CASE(1)
+              ALLOCATE(NEW_LIST_DP1(newSize),STAT=ERR)
+              IF(ERR/=0) CALL FlagError("Could not allocate new list.",ERR,ERROR,*999)
+              NEW_LIST_DP1(1:LIST%NUMBER_IN_LIST)=LIST%LIST_DP(1:LIST%NUMBER_IN_LIST)
+              CALL MOVE_ALLOC(NEW_LIST_DP1,LIST%LIST_DP)
+            CASE(2:) 
+              ALLOCATE(NEW_LIST_DP2(LIST%DATA_DIMENSION,newSize),STAT=ERR)
+              IF(ERR/=0) CALL FlagError("Could not allocate new list.",ERR,ERROR,*999)
+              NEW_LIST_DP2(:,1:LIST%NUMBER_IN_LIST)=LIST%LIST_DP2(:,1:LIST%NUMBER_IN_LIST)
+              CALL MOVE_ALLOC(NEW_LIST_DP2,LIST%LIST_DP2)
+            CASE DEFAULT
+              CALL FlagError("Invalid data dimension.",ERR,ERROR,*999)
+            END SELECT
+          CASE DEFAULT
+            CALL FlagError("Invalid data type.",ERR,ERROR,*999)
+          END SELECT
+            
+          LIST%SIZE=newSize
+
+        ENDIF
+      ELSE
+        CALL FlagError("The list has not been finished.",ERR,ERROR,*999)
+      ENDIF
+    ELSE
+      CALL FlagError("List is not associated.",ERR,ERROR,*999)
+    ENDIF
+    
+    EXITS("List_Resize")
+    RETURN
+999 IF(ALLOCATED(NEW_LIST_INT1)) DEALLOCATE(NEW_LIST_INT1)
+    IF(ALLOCATED(NEW_LIST_INT2)) DEALLOCATE(NEW_LIST_INT2)
+    IF(ALLOCATED(NEW_LIST_SP1)) DEALLOCATE(NEW_LIST_SP1)
+    IF(ALLOCATED(NEW_LIST_SP2)) DEALLOCATE(NEW_LIST_SP2)
+    IF(ALLOCATED(NEW_LIST_DP1)) DEALLOCATE(NEW_LIST_DP1)
+    IF(ALLOCATED(NEW_LIST_DP2)) DEALLOCATE(NEW_LIST_DP2)
+    ERRORSEXITS("List_Resize",ERR,ERROR)
+    RETURN 1
+    
+  END SUBROUTINE List_Resize
+
+
   !
   !================================================================================================================================
   !
@@ -979,6 +1085,7 @@ SUBROUTINE LIST_ITEM_ADD_INTG2(LIST,ITEM,ERR,ERROR,*)
           IF(LIST%DATA_DIMENSION==SIZE(ITEM,1)) THEN
             IF(LIST%NUMBER_IN_LIST==LIST%SIZE) THEN
               !Reallocate
+              ! DOUBLE the existing list size if list is full
               NEW_SIZE=MAX(2*LIST%NUMBER_IN_LIST,1)
               ALLOCATE(NEW_LIST(LIST%DATA_DIMENSION,NEW_SIZE),STAT=ERR)
               IF(ERR/=0) CALL FlagError("Could not allocate new list.",ERR,ERROR,*999)
diff --git a/src/mesh_routines.F90 b/src/mesh_routines.F90
index e5d467f7..96ed5ef4 100644
--- a/src/mesh_routines.F90
+++ b/src/mesh_routines.F90
@@ -855,17 +855,25 @@ SUBROUTINE DECOMPOSITION_ELEMENT_DOMAIN_GET(DECOMPOSITION,USER_ELEMENT_NUMBER,DO
                 CALL TREE_NODE_VALUE_GET(MESH_ELEMENTS%ELEMENTS_TREE,TREE_NODE,GLOBAL_ELEMENT_NUMBER,ERR,ERROR,*999)
                 IF(GLOBAL_ELEMENT_NUMBER>0.AND.GLOBAL_ELEMENT_NUMBER<=MESH_TOPOLOGY%ELEMENTS%NUMBER_OF_ELEMENTS) THEN
                   DOMAIN_NUMBER=DECOMPOSITION%ELEMENT_DOMAIN(GLOBAL_ELEMENT_NUMBER)
+                
+                    ! adding an output here!!!
+                    !IF (DOMAIN_NUMBER==0) THEN
+                    !  CALL TREE_OUTPUT(1,MESH_ELEMENTS%ELEMENTS_TREE,ERR,ERROR,*999)
+                    !END IF
+
                 ELSE
                   LOCAL_ERROR="Global element number found "//TRIM(NUMBER_TO_VSTRING(GLOBAL_ELEMENT_NUMBER,"*",ERR,ERROR))// &
                     & " is invalid. The limits are 1 to "// &
                     & TRIM(NUMBER_TO_VSTRING(MESH_TOPOLOGY%ELEMENTS%NUMBER_OF_ELEMENTS,"*",ERR,ERROR))//"."
-                  CALL FlagError(LOCAL_ERROR,ERR,ERROR,*999)
+                  !CALL FlagError(LOCAL_ERROR,ERR,ERROR,*999)
                 ENDIF
               ELSE
                 LOCAL_ERROR="Decomposition mesh element corresponding to user element number "// &
                   & TRIM(NumberToVString(USER_ELEMENT_NUMBER,"*",err,error))//" is not found."
-                CALL FlagError(LOCAL_ERROR,ERR,ERROR,*999)
-              ENDIF
+                  ! CALL FlagError(LOCAL_ERROR,ERR,ERROR,*999)
+                  ! restore this flag!!!! 
+                  ERR = 1000
+             ENDIF
             ELSE
               CALL FlagError("Decomposition mesh elements are not associated.",ERR,ERROR,*999)
             ENDIF
@@ -1657,6 +1665,13 @@ SUBROUTINE DecompositionTopology_ElementCheckExists(decompositionTopology,userEl
       CALL Tree_NodeValueGet(decompositionElements%ELEMENTS_TREE,treeNode,localElementNumber,err,error,*999)
       elementExists=.TRUE.
       ghostElement=localElementNumber>decompositionElements%NUMBER_OF_ELEMENTS
+
+                    ! how to call the topology tree?...
+                    !IF (DOMAIN_NUMBER==0) THEN
+                    ! PRINT *, "Topology tree "
+                    ! CALL TREE_OUTPUT(1,decompositionElements%ELEMENTS_TREE,ERR,ERROR,*999)
+                    !END IF
+
     ENDIF
     
     EXITS("DecompositionTopology_ElementCheckExists")
diff --git a/src/opencmiss_iron.F90 b/src/opencmiss_iron.F90
index ed0b281f..ac674852 100644
--- a/src/opencmiss_iron.F90
+++ b/src/opencmiss_iron.F90
@@ -29250,6 +29250,8 @@ SUBROUTINE cmfe_Field_CreateFinishObj(field,err)
     CALL TAU_STATIC_PHASE_STOP('field Create')
 #endif
 
+! another comment
+
     EXITS("cmfe_Field_CreateFinishObj")
     RETURN
 999 ERRORSEXITS("cmfe_Field_CreateFinishObj",err,error)
diff --git a/src/solver_mapping_routines.F90 b/src/solver_mapping_routines.F90
index 4e680eae..6ccf87de 100644
--- a/src/solver_mapping_routines.F90
+++ b/src/solver_mapping_routines.F90
@@ -314,7 +314,7 @@ SUBROUTINE SOLVER_MAPPING_CALCULATE(SOLVER_MAPPING,ERR,ERROR,*)
           !
           !--- Row mappings ---
           !
-          ! 2. Determine the number of rows in the solver matrix. Do this the by setting up a list of rows for each rank.
+          ! 2. Determine the number of rows in the solver matrix. Do this by setting up a list of rows for each rank.
           !    We can then later arrange the rows in rank order by looping over the ranks in the list and then the rows
           !    for each rank.
           !
@@ -395,9 +395,9 @@ SUBROUTINE SOLVER_MAPPING_CALCULATE(SOLVER_MAPPING,ERR,ERROR,*)
               ENDIF
             ENDDO !interface_condition_idx
           ENDDO !rank
-          !Calculate the number of local and global rows. Do this by looking at the boundary conditions for field variables
-          !involved in the row. If all the variables are set as a fixed boundary condition then do not include the row. If
-          !any variable is not fixed then include the row.
+          !Calculate the NUMBER of local and global rows. Do this by looking at the boundary conditions for field variables
+          !involved in the row. ONLY If ALL the variables are set as a fixed boundary condition then do NOT include the row. If
+          !any variable is not fixed then INCLUDE the row.
           equations_idx=0
           DO equations_set_idx=1,SOLVER_MAPPING%NUMBER_OF_EQUATIONS_SETS
             equations_idx=equations_idx+1
@@ -426,6 +426,7 @@ SUBROUTINE SOLVER_MAPPING_CALCULATE(SOLVER_MAPPING,ERR,ERROR,*)
                         IF(ROW_DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(global_row)%LOCAL_TYPE(rank_idx)/=DOMAIN_LOCAL_GHOST) THEN
                           ROW_RANK=ROW_DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(global_row)%DOMAIN_NUMBER(rank_idx)
                           local_row=ROW_DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(global_row)%LOCAL_NUMBER(rank_idx)
+                        ! EXCLUDE rows that are only ghosts!!!
                           EXIT
                         ENDIF
                       ENDDO !rank_idx
@@ -440,7 +441,9 @@ SUBROUTINE SOLVER_MAPPING_CALCULATE(SOLVER_MAPPING,ERR,ERROR,*)
                           IF(ASSOCIATED(BOUNDARY_CONDITIONS_VARIABLE)) THEN
                             !This is wrong as we only have the mappings for the local rank not the global ranks.
                             !For now assume 1-1 mapping between rows and dofs.
-                            global_dof=global_row
+                            ! In progress: Consider eigentransformations for constrained equations.
+                            global_dof=global_row ! assumption!!!
+
                             INCLUDE_ROW=INCLUDE_ROW.AND.(BOUNDARY_CONDITIONS_VARIABLE%DOF_TYPES(global_dof)== &
                               & BOUNDARY_CONDITION_DOF_FREE)
                             CONSTRAINED_DOF=CONSTRAINED_DOF.OR.(BOUNDARY_CONDITIONS_VARIABLE%DOF_TYPES(global_dof)== &
@@ -562,7 +565,10 @@ SUBROUTINE SOLVER_MAPPING_CALCULATE(SOLVER_MAPPING,ERR,ERROR,*)
               CALL FlagError("Equations set is not associated.",ERR,ERROR,*999)
             ENDIF
           ENDDO !equations set idx
+
+          ! DO THE SAME FOR INTERFACE BLOCKS!!!
           !Now add in rows from any interface matrices
+          !WHICH are transpose, THEN look at the columns!!!
           DO interface_condition_idx=1,SOLVER_MAPPING%NUMBER_OF_INTERFACE_CONDITIONS
             equations_idx=equations_idx+1
             INTERFACE_CONDITION=>SOLVER_MAPPING%INTERFACE_CONDITIONS(interface_condition_idx)%PTR
@@ -745,7 +751,7 @@ SUBROUTINE SOLVER_MAPPING_CALCULATE(SOLVER_MAPPING,ERR,ERROR,*)
               CALL SolverMapping_InterfToSolMatMapsIMInitialise(SOLVER_MAPPING%INTERFACE_CONDITION_TO_SOLVER_MAP( &
                 & interface_condition_idx)%INTERFACE_TO_SOLVER_MATRIX_MAPS_IM(interface_matrix_idx),ERR,ERROR,*999)
                       
-              !Allocate the interfafce row to solver row maps
+              !Allocate the interface row to solver row maps
               ALLOCATE(SOLVER_MAPPING%INTERFACE_CONDITION_TO_SOLVER_MAP(interface_condition_idx)% &
                 & INTERFACE_TO_SOLVER_MATRIX_MAPS_IM(interface_matrix_idx)%INTERFACE_ROW_TO_SOLVER_ROWS_MAP( &
                 & INTERFACE_MAPPING%INTERFACE_MATRIX_ROWS_TO_VAR_MAPS(interface_matrix_idx)%TOTAL_NUMBER_OF_ROWS),STAT=ERR)
@@ -760,7 +766,7 @@ SUBROUTINE SOLVER_MAPPING_CALCULATE(SOLVER_MAPPING,ERR,ERROR,*)
               
             ENDDO !interface_matrix_idx
 
-            !Allocate the interface column to solver row maps
+            !Allocate the interface column(!!!!transpose!!!!!!) to solver row maps
             ALLOCATE(SOLVER_MAPPING%INTERFACE_CONDITION_TO_SOLVER_MAP(interface_condition_idx)% &
               & INTERFACE_COLUMN_TO_SOLVER_ROWS_MAPS(INTERFACE_MAPPING%TOTAL_NUMBER_OF_COLUMNS),STAT=ERR)
             IF(ERR/=0)  &
@@ -841,6 +847,8 @@ SUBROUTINE SOLVER_MAPPING_CALCULATE(SOLVER_MAPPING,ERR,ERROR,*)
                   !Set up the row domain mappings.
                   !There are no ghosted rows for the solver matrices so there is only one domain for the global to local map.
                   !Initialise
+                  ! THESE global to local will have to be replaced with MPI communication (reduction)!!!!!!!!
+                  ! Knowledge of each rank about the rows in other ranks!
                   CALL DOMAIN_MAPPINGS_MAPPING_GLOBAL_INITIALISE(ROW_DOMAIN_MAPPING%GLOBAL_TO_LOCAL_MAP( &
                     & NUMBER_OF_GLOBAL_SOLVER_ROWS),ERR,ERROR,*999)
                   !Allocate the global to local map arrays
diff --git a/src/trees.F90 b/src/trees.F90
index 9e994023..165c6d7e 100644
--- a/src/trees.F90
+++ b/src/trees.F90
@@ -100,7 +100,7 @@ MODULE TREES
     INTEGER(INTG) :: INSERT_TYPE !<The insert type for duplicate keys for the tree
     INTEGER(INTG) :: NUMBER_IN_TREE !<The number of items currently in the tree
     TYPE(TREE_NODE_TYPE), POINTER :: ROOT !<The pointer to the root of the tree
-    TYPE(TREE_NODE_TYPE), POINTER :: NIL !<The pointer to the nil of the tree
+    TYPE(TREE_NODE_TYPE), POINTER :: NIL !<The pointer to the nil (= leaves) of the tree
   END TYPE TREE_TYPE
 
   !Module variables
@@ -287,13 +287,13 @@ END SUBROUTINE TREE_DESTROY
   !================================================================================================================================
   !
 
-  !>Detaches the tree values and returns them as a pointer to the an array
+  !>Detaches the tree values and returns them as a pointer to an array
   SUBROUTINE TREE_DETACH(TREE,NUMBER_IN_TREE,TREE_VALUES,ERR,ERROR,*)
 
     !Argument Variables
     TYPE(TREE_TYPE), POINTER :: TREE !<A pointer to the tree to detach
     INTEGER(INTG), INTENT(OUT) :: NUMBER_IN_TREE !<On exit, the number in the array that has been detached
-    INTEGER(INTG), POINTER :: TREE_VALUES(:) !<On exit, a pointer to the dettached tree values. Must not be associated on entry.
+    INTEGER(INTG), POINTER :: TREE_VALUES(:) !<On exit, a pointer to the detached tree values. Must not be associated on entry.
     INTEGER(INTG), INTENT(OUT) :: ERR !<The error code
     TYPE(VARYING_STRING), INTENT(OUT) :: ERROR !<The error string.
     !Local Variables
@@ -330,7 +330,7 @@ END SUBROUTINE TREE_DETACH
   !================================================================================================================================
   !
 
-  !>Detaches the tree values and returns them as a pointer to the an array and then destroys the tree
+  !>Detaches the tree values and returns them as a pointer to an array and then destroys the tree
   SUBROUTINE TREE_DETACH_AND_DESTROY(TREE,NUMBER_IN_TREE,TREE_VALUES,ERR,ERROR,*)
 
     !Argument Variables
@@ -374,7 +374,7 @@ END SUBROUTINE TREE_DETACH_AND_DESTROY
   !================================================================================================================================
   !
 
-  !>Detaches the tree values in order from the specified tree node and adds them to the tree values array
+  !>Detaches the tree values in order from the specified tree node (eg root) and adds them to the tree values array
   RECURSIVE SUBROUTINE TREE_DETACH_IN_ORDER(TREE,X,COUNT,TREE_VALUES,ERR,ERROR,*)
 
     !Argument Variables
diff --git a/src/types.F90 b/src/types.F90
index f40a2c50..52179e49 100644
--- a/src/types.F90
+++ b/src/types.F90
@@ -97,7 +97,8 @@ MODULE Types
   TYPE INTEGER_CINT_ALLOC_TYPE
     INTEGER(C_INT), ALLOCATABLE :: ARRAY(:)
   END TYPE INTEGER_CINT_ALLOC_TYPE
-  
+
+
   !
   !================================================================================================================================
   !
@@ -130,6 +131,32 @@ MODULE Types
     INTEGER(C_INT), ALLOCATABLE :: LIST_C_INT(:) !<The integer data (dimension = 1) for integer lists. 
     INTEGER(C_INT), ALLOCATABLE :: LIST_C_INT2(:,:) !<The integer data (dimension > 1) for integer lists. 
   END TYPE LIST_TYPE
+
+
+  !
+  !================================================================================================================================
+  !
+  ! Hash Table type
+  !
+  TYPE HashTableType
+    !PRIVATE ! derived-type definition can be accessed outside (public), but not the components (below)
+
+    LOGICAL :: hashTableFinished !<Is .TRUE. if the table has finished being created, .FALSE. if not.
+
+    INTEGER(INTG), ALLOCATABLE :: vecTKey(:)   ! The hash vector of keys
+    INTEGER(INTG), ALLOCATABLE :: vecTVal(:)   ! The hash vector of indices of values
+
+    INTEGER(INTG), ALLOCATABLE :: vecSKey(:)   ! The original vector of keys
+ 
+    ! List of values
+    TYPE(List_Type), POINTER :: listSVal
+    ! Array of list of values to allow for multiple data
+    TYPE(LIST_PTR_TYPE), ALLOCATABLE :: arrayOfListSVal(:)
+
+
+    INTEGER(INTG) :: n !<The number of items currently in the table (number of keys)
+    INTEGER(INTG) :: p !<The prime number required from the table algorithm
+  END TYPE HashTableType
     
   !
   !================================================================================================================================
@@ -847,7 +874,10 @@ MODULE Types
   TYPE DistributedMatrixCMISSType
     TYPE(DistributedMatrixType), POINTER :: distributedMatrix !<A pointer to the distributed matrix
     INTEGER(INTG) :: baseTagNumber !<The base number for the MPI tag numbers that will be used to communicate the distributed matrix data amongst the domains. The base tag number can be thought of as the identification number for the distributed matrix object.
-    TYPE(MATRIX_TYPE), POINTER :: matrix !<A pointer to the matrix to store the rows corresponding to this domain.
+    TYPE(MATRIX_TYPE), POINTER   :: matrix !<A pointer to the matrix to store the rows corresponding to this domain.
+    
+    TYPE(HashTableType), POINTER :: columnHashTable !< The hash table for ltg columns
+
   END TYPE DistributedMatrixCMISSType
 
   !>Contains information for a PETSc distributed matrix
@@ -1998,8 +2028,11 @@ MODULE Types
     TYPE(BOUNDARY_CONDITIONS_TYPE), POINTER :: BOUNDARY_CONDITIONS !<A pointer to the boundary conditions for this boundary conditions variable
     INTEGER(INTG) :: VARIABLE_TYPE !<The type of variable for this variable boundary conditions
     TYPE(FIELD_VARIABLE_TYPE), POINTER :: VARIABLE !<A pointer to the field variable for this boundary condition variable
+
+! These two types below have been confused!!!
     INTEGER(INTG), ALLOCATABLE :: DOF_TYPES(:) !<DOF_TYPES(dof_idx). The general boundary condition type (eg. fixed or free) of the dof_idx'th dof in the dependent field variable. \see OPENCMISS_BoundaryConditionsTypes,OPENCMISS
     INTEGER(INTG), ALLOCATABLE :: CONDITION_TYPES(:) !<CONDITION_TYPES(dof_idx). The specific boundary condition type (eg. incremented pressure) of the dof_idx'th dof of the dependent field variable, which might be specific to an equation set. The solver routines should not need to use this array, and should only need the DOF_TYPES array. \see OPENCMISS_BoundaryConditionsDOFTypes,OPENCMISS
+    
     TYPE(BOUNDARY_CONDITIONS_DIRICHLET_TYPE), POINTER :: DIRICHLET_BOUNDARY_CONDITIONS  !<A pointer to the dirichlet boundary condition type for this boundary condition variable
     INTEGER(INTG) :: NUMBER_OF_DIRICHLET_CONDITIONS !<Stores the number of dirichlet conditions associated with this variable
     TYPE(BoundaryConditionsNeumannType), POINTER :: neumannBoundaryConditions
@@ -2045,6 +2078,8 @@ MODULE Types
   !>Contains information used to integrate Neumann boundary conditions
   TYPE BoundaryConditionsNeumannType
     INTEGER(INTG), ALLOCATABLE :: setDofs(:) !<setDofs(neumann_idx): the global dof for the neumann_idx'th Neumann condition
+
+    ! This matrix should be local!!!??
     TYPE(DistributedMatrixType), POINTER :: integrationMatrix !<The N matrix that multiples the point values vector q to give the integrated values f. Number of rows equals number of local dofs, and number of columns equals number of set point DOFs.
     TYPE(DistributedVectorType), POINTER :: pointValues !<The vector of set point values q
     TYPE(DOMAIN_MAPPING_TYPE), POINTER :: pointDofMapping !<The domain mapping for DOFs with Neumann point conditions set.