anchor.diff:      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
  +       num_elems_new = 2*num_elems
       elseif(mesh_type.eq.'terminal')then
  -      num_elems_new = 2*num_elems + num_units
  +       num_elems_new = 2*num_elems + num_units
       elseif(mesh_type.eq.'ladder')then
  -        print *, "NOT YET IMPLEMENTED"
  -        call exit(0)
  +       print *, "NOT YET IMPLEMENTED"
  +       call exit(0)
       endif
       call reallocate_node_elem_arrays(num_elems_new,num_nodes_new)
       noelem0=0
  
positive.code:     if(allocated(elem_direction)) deallocate(elem_direction)
      allocate(elem_direction(3,num_elems))
      if(model_type.eq.'gas_mix')then
        if(allocated(expansile)) deallocate(expansile)
        allocate(expansile(num_elems))
      endif
  
  !!! initialise element arrays
      elems=0
      elem_nodes=0
      elem_symmetry = 1
      elem_field = 0.0_dp
      if(model_type.eq.'gas_mix')expansile = .false.
  
      ne=0
  
      read_an_element : do
         !.......read element number
         read(unit=10, fmt="(a)", iostat=ierror) ctemp1
  
negative.code:             elem_field(ne_length,ne)
         elem_field(ne_a_A,ne) = 1.0_dp ! set default for ratio a/A
      enddo
      
      call element_connectivity_1d
      call evaluate_ordering
      
      call enter_exit(sub_name,2)
      
    end subroutine define_mesh_geometry_test
  
  !!!#############################################################################
    
    subroutine define_node_geometry(NODEFILE)
      !*define_node_geometry:* Reads in an ipnode file to define a tree geometry
      !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_DEFINE_NODE_GEOMETRY" :: DEFINE_NODE_GEOMETRY
      
      character(len=MAX_FILENAME_LEN), intent(in) :: NODEFILE !Input nodefile
      !     Local Variables
      integer :: i,ierror,np,np_global,num_nodes_temp,num_versions,nv,NJT=0
      real(dp) :: point
      logical :: overwrite = .false. ! initialised
      character(len=300) :: ctemp1,readfile
      character(len=60) :: sub_name
      
      ! --------------------------------------------------------------------------
      
      sub_name = 'define_node_geometry'
      call enter_exit(sub_name,1)
      
      readfile = trim(nodefile)//'.ipnode'
      open(10, file=NODEFILE, status='old')
      
      if(num_nodes.gt.0) overwrite = .true.
      
      !.....read in the total number of nodes. read each line until one is found
  