This routine contains a loop over all edges. For each edge i which is to be divided into two sub-edges the following actions are taken:

- determine the coordinates of the new node on the edge and store them in a consecutively arranged field conewnodes(1..3,*) with the total number of new nodes stored in the scalar nnewnodes. Let us assume that edge i is the j-th edge to be subdivided, then the coordinates are stored in conewnodes(1..3,j)
- determine a base element to which the edge belongs (an arbitrary element of the shell of the edge) and store the number in ibasenewnodes(j)
- store the edge number on which the nodes lies (i.e. here i) in field iedgnewnodes(*), i.e. iedgnewnodes(j)=i.
- determine the value of the h-field in the new node by interpolation within the unrefined mesh and store it in hnewnodes(j)

After leaving newnodes a random contribution is added to the coordinates of each new node while moving them towards the center of gravity of their base element (=perturbation of the nodal position from a surface position to a subsurface position). This ensures that each new node lies within a tetrahedral element and not on its faces or edges. This facilitates the insertion of the new nodes in the existing mesh and is particularly important for the newly created external surface nodes. Indeed, the insertion procedure explained in [24] works only for nodes not lying on the external surface. The insertion of the new nodes is done one by one after reordering them in an aleatoric way.