anchor.diff:      call enter_exit(sub_name,2)
   
     end subroutine reallocate_node_elem_arrays
  +
   !
   !###################################################################################
   !
  -!>define_mesh_geometry_test
  -  subroutine define_mesh_geometry_test
  -    use arrays,only: dp,nodes,node_field,node_xyz,num_nodes,&
  -         elem_direction,elem_field,elems,elem_cnct,elem_nodes,&
  -         elem_ordrs,elem_symmetry,elems_at_node,elem_units_below,&
  -         expansile,node_xyz,num_elems,num_nodes
  -    use indices
  -    implicit none
  -
  -    !     Local Variables
  -    integer :: j,ne,np,np1,np2
  -
  -    num_elems = 400
  -    num_nodes = num_elems + 1
  -
  -!!! allocate memory
  -    if(.not.allocated(nodes)) allocate (nodes(num_nodes))
  -    if(.not.allocated(node_xyz)) allocate (node_xyz(3,num_nodes))
  -    if(.not.allocated(node_field)) allocate (node_field(num_nj,num_nodes))
  -    if(.not.allocated(elems)) allocate(elems(num_elems))
  -    if(.not.allocated(elem_cnct)) allocate(elem_cnct(-1:1,0:2,0:num_elems))
  -    if(.not.allocated(elem_nodes)) allocate(elem_nodes(2,num_elems))
  -    if(.not.allocated(elem_ordrs)) allocate(elem_ordrs(num_ord,num_elems))
  -    if(.not.allocated(elem_symmetry)) allocate(elem_symmetry(num_elems))
  -    if(.not.allocated(elem_units_below)) allocate(elem_units_below(num_elems))
  -    if(.not.allocated(elems_at_node)) allocate(elems_at_node(num_nodes,0:3))
  -    if(.not.allocated(elem_field)) allocate(elem_field(num_ne,num_elems))
  -    if(.not.allocated(elem_direction)) allocate(elem_direction(3,num_elems))
  -    if(.not.allocated(expansile)) allocate(expansile(num_elems))
  -
  -!!! initialise array values
  -    nodes = 0 !initialise node index values
  -    node_xyz = 0.0_dp !initialise
  -    node_field = 0.0_dp !initialise
  -    elems=0
  -    elem_nodes=0
  -    elem_symmetry = 1
  -    elem_field = 0.0_dp
  -    expansile = .false.
  -
  -!!! set up node arrays
  -    nodes(1) = 1
  -    do np=2,101
  -       nodes(np) = np
  -       node_xyz(3,np) = node_xyz(3,np-1) - 1.0_dp
  -    enddo
  -
  -    np=102
  -    node_xyz(1,np) = node_xyz(1,np-1) - 0.5_dp
  -    do np=102,151
  -       nodes(np) = np
  -       node_xyz(3,np) = node_xyz(3,np-1) - 0.5_dp
  -       node_xyz(1,np) = node_xyz(1,np-1) - 0.5_dp
  -    enddo
  -
  -    np=152
  -    nodes(np) = np
  -    node_xyz(1,np) = node_xyz(1,101) + 0.5_dp
  -    node_xyz(3,np) = node_xyz(3,101) - 0.5_dp
  -    do np=153,201
  -       nodes(np) = np
  -       node_xyz(3,np) = node_xyz(3,np-1) - 0.5_dp
  -       node_xyz(1,np) = node_xyz(1,np-1) + 0.5_dp
  -    enddo
  -
  -    np=202
  -    nodes(np) = np
  -    node_xyz(1,np) = node_xyz(1,151) - 0.5_dp
  -    node_xyz(3,np) = node_xyz(3,151) - 0.5_dp
  -    do np=203,251
  -       nodes(np) = np
  -       node_xyz(3,np) = node_xyz(3,np-1) - 0.5_dp
  -       node_xyz(1,np) = node_xyz(1,np-1) - 0.5_dp
  -    enddo
  -
  -    np=252
  -    nodes(np) = np
  -    node_xyz(1,np) = node_xyz(1,151) + 0.5_dp
  -    node_xyz(3,np) = node_xyz(3,151) - 0.5_dp
  -    do np=253,301
  -       nodes(np) = np
  -       node_xyz(3,np) = node_xyz(3,np-1) - 0.5_dp
  -       node_xyz(1,np) = node_xyz(1,np-1) + 0.5_dp
  -    enddo
  -
  -    np=302
  -    nodes(np) = np
  -    node_xyz(1,np) = node_xyz(1,201) - 0.5_dp
  -    node_xyz(3,np) = node_xyz(3,201) - 0.5_dp
  -    do np=303,351
  -       nodes(np) = np
  -       node_xyz(3,np) = node_xyz(3,np-1) - 0.5_dp
  -       node_xyz(1,np) = node_xyz(1,np-1) - 0.5_dp
  -    enddo
  -
  -    np=352
  -    nodes(np) = np
  -    node_xyz(1,np) = node_xyz(1,201) + 0.5_dp
  -    node_xyz(3,np) = node_xyz(3,201) - 0.5_dp
  -    do np=353,401
  -       nodes(np) = np
  -       node_xyz(3,np) = node_xyz(3,np-1) - 0.5_dp
  -       node_xyz(1,np) = node_xyz(1,np-1) + 0.5_dp
  -    enddo
  -
  -!!! set up element arrays
  -    do ne=1,num_elems
  -       elems(ne) = ne
  -       elem_nodes(1,ne) = ne
  -       elem_nodes(2,ne) = ne+1
  -    enddo
  -
  -    elem_nodes(1,151) = 101
  -    elem_nodes(1,201) = 151
  -    elem_nodes(1,251) = 151
  -    elem_nodes(1,301) = 201
  -    elem_nodes(1,351) = 201
  -
  -    elem_field(ne_radius,1:100) = 10.0_dp
  -    elem_field(ne_radius,101:200) = 5.0_dp
  -    elem_field(ne_radius,201:400) = dsqrt(elem_field(ne_radius,101)**2/2.0_dp)
  -
  -    ! calculate the element lengths and directions
  -    do ne=1,num_elems
  -       np1=elem_nodes(1,ne)
  -       np2=elem_nodes(2,ne)
  -       elem_field(ne_length,ne) = DSQRT((node_xyz(1,np2) - &
  -            node_xyz(1,np1))**2 + (node_xyz(2,np2) - &
  -            node_xyz(2,np1))**2 + (node_xyz(3,np2) - &
  -            node_xyz(3,np1))**2)
  -       do j=1,3
  -          elem_direction(j,ne) = (node_xyz(j,np2) - &
  -               node_xyz(j,np1))/elem_field(ne_length,ne)
  -       enddo !j
  -       ! element volume
  -       elem_field(ne_vol,ne) = PI * elem_field(ne_radius,ne)**2 * &
  -            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
  -
  -  end subroutine define_mesh_geometry_test
  -
  -!
  -!######################################################################################################
  +!>get_final_real
  +  subroutine get_final_real(string,rtemp)
  +    use arrays,only: dp
  +    implicit none
  +    character, intent(in) :: string*(132)
  +    integer :: ibeg,iend
  +    real(dp), intent(out) :: rtemp
  +    real(dp) :: rsign
  +    character :: sub_string*(40)
  +
  +    iend=len(string) !get the length of the string
  +    ibeg=index(string,":")+1 !get location of real in string
  +    sub_string = adjustl(string(ibeg:iend)) ! get the characters beyond :
  +    iend=len(sub_string) !get the length of the sub-string
  +    if(sub_string(1:1).eq.'-')then !check whether negative
  +       rsign=-1.0_dp
  +       ibeg=2
  +    else
  +       rsign=1.0_dp
  +       ibeg=1
  +    endif
  +    read (sub_string(ibeg:iend), '(D25.17)' ) rtemp !get real value
  +    rtemp=rtemp*rsign !apply sign to number
  +
  +  end subroutine get_final_real
  +
  +!
  +!###################################################################################
  +!
  +!*get_final_string*
  +  subroutine get_final_string(string,rtemp)
  +    use arrays,only: dp
  +    implicit none
  +    character, intent(in) :: string*(132)
  +    integer :: ibeg,iend
  +    real(dp), intent(out) :: rtemp
  +    real(dp) :: rsign
  +    character :: sub_string*(40)
  +
  +    iend=len(string) !get the length of the string
  +    ibeg=index(string,":")+1 !get location of real in string
  +    sub_string = adjustl(string(ibeg:iend)) ! get the characters beyond :
  +    iend=len(sub_string) !get the length of the sub-string
  +    if(sub_string(1:1).eq.'-')then !check whether negative
  +       rsign=-1.0_dp
  +       ibeg=2
  +    else
  +       rsign=1.0_dp
  +       ibeg=1
  +    endif
  +    read (sub_string(ibeg:iend), '(D25.17)' ) rtemp !get real value
  +    rtemp=rtemp*rsign !apply sign to number
  +
  +  end subroutine get_final_string
  +
  +!
  +!###################################################################################
  +!
  +!*get_local_node*
  +  subroutine get_local_node(np_global,np_local)
  +    use arrays,only: nodes,num_nodes
  +    implicit none
  +
  +    integer,intent(in) :: np_global
  +    integer,intent(out) :: np_local
  +
  +    integer :: np
  +    logical :: found
  +
  +    np=1
  +    found=.false.
  +    do while (.not.found)
  +       if(nodes(np).eq.np_global)then
  +          found=.true.
  +       elseif(np.gt.num_nodes)then
  +          found = .true.
  +          write(*,'('' Global node '',I6,'' not in node list'')') np_global
  +          read(*,*)
  +       else
  +          np=np+1
  +       endif
  +    enddo
  +
  +    np_local = np
  +    return
  +
  +  end subroutine get_local_node
  +
  +!
  +!###################################################################################
  +!
  +!*get_final_integer*
  +  subroutine get_final_integer(string,num)
  +    implicit none
  +    character,intent(in) :: string*(132)
  +    integer,intent(out) :: num
  +    integer :: ibeg,iend,nsign,ntemp
  +    character :: sub_string*(40)
  +
  +    iend=len(string) !get the length of the string
  +    ibeg=index(string,":")+1 !get location of integer in string, follows ":"
  +    sub_string = adjustl(string(ibeg:iend)) ! get the characters beyond ":"
  +    iend=len(sub_string) !length of the sub-string
  +    if(sub_string(1:1).eq.'-')then !check for negative sign
  +       nsign=-1
  +       ibeg=2
  +    else
  +       nsign=1
  +       ibeg=1
  +    endif
  +    read (sub_string(ibeg:iend), '(i10)' ) ntemp !get integer values
  +    ntemp=ntemp*nsign !apply sign to number
  +
  +    num=ntemp !return the integer value
  +
  +  end subroutine get_final_integer
  +
  +!
  +!###########################################################################################
   !
   end module geometry
   
  
positive.code:     call enter_exit(sub_name,2)
  
    END subroutine define_rad_from_file
  
  !
  !##################################################################################
  !
  !>define_rad_from_geom
    subroutine define_rad_from_geom(ORDR_SYSTEM,CONTROL_PARAM,START_FROM,START_RAD,GROUP_TYPE,GROUP_OPTIONS)
      use arrays,only: dp,num_elems,elem_field,elem_ordrs,maxgen,elem_cnct
      use indices
  
negative.code:    integer :: ne_min,ne_max,nindex,ne,n_max_ord,n,ne_start,&
        inlet_count
     real(dp) :: radius
     character(len=60) :: sub_name
  
     sub_name = 'define_rad_from_geom'
  