anchor.diff:      elem_direction(1:3,1:num_elems)=rnodes_temp(1:3,1:num_elems)
       deallocate(rnodes_temp)
       elem_direction(1:3,num_elems+1:num_elems_new) = 0.0_dp
  -
  +    
       if(allocated(node_field).and.num_nj.gt.0)then
  -      allocate(rnodes_temp(num_nj,num_nodes))
  -      rnodes_temp=node_field
  -      deallocate(node_field)
  -      allocate(node_field(num_nj,num_nodes_new))
  -      node_field(1:num_nj,1:num_nodes)=rnodes_temp(1:num_nj,1:num_nodes)
  -      deallocate(rnodes_temp)
  -      node_field(1:num_nj,num_nodes+1:num_nodes_new)=0.0_dp
  +       allocate(rnodes_temp(num_nj,num_nodes))
  +       rnodes_temp=node_field
  +       deallocate(node_field)
  +       allocate(node_field(num_nj,num_nodes_new))
  +       node_field(1:num_nj,1:num_nodes)=rnodes_temp(1:num_nj,1:num_nodes)
  +       deallocate(rnodes_temp)
  +       node_field(1:num_nj,num_nodes+1:num_nodes_new)=0.0_dp
       endif
  -
  +    
       allocate(nodelem_temp(num_elems))
       nodelem_temp = elem_symmetry ! copy to temporary array
       deallocate(elem_symmetry) !deallocate initially allocated memory
  
positive.code:     num_nodes_new = 2*num_nodes
      ! the number of elems after adding mesh will be:
      if(mesh_type.eq.'basic')then
        num_elems_new = 2*num_elems
      elseif(mesh_type.eq.'terminal')then
        num_elems_new = 2*num_elems + num_units
      elseif(mesh_type.eq.'ladder')then
          print *, "NOT YET IMPLEMENTED"
          call exit(0)
      endif
      call reallocate_node_elem_arrays(num_elems_new,num_nodes_new)
      noelem0=0
  
negative.code:           exit read_versions !exit the named do loop
         endif
      end do read_versions
      
      np = 0
      !.....read the coordinate, derivative, and version information for each node.
      read_a_node : do !define a do loop name
  