anchor.diff:  
     !Interfaces
     private
  -  public add_mesh,define_rad_from_file, define_rad_from_geom,&
  -       define_node_geometry,&
  -       define_1d_elements,element_connectivity_1d,&
  -       evaluate_ordering,volume_of_mesh,set_initial_volume,&
  -       append_units,define_mesh_geometry_test
  +  public add_mesh
  +  public append_units
  +  public define_1d_elements
  +  public define_mesh_geometry_test
  +  public define_node_geometry
  +  public define_rad_from_file
  +  public define_rad_from_geom
  +  public element_connectivity_1d
  +  public evaluate_ordering
  +  public set_initial_volume
  +  public volume_of_mesh
   
   contains
  +!
  +!###################################################################################
  +!
  +!*add_mesh:* Reads in an ipmesh file and adds this mesh to the terminal branches of an existing tree geometry
  +  subroutine add_mesh(AIRWAY_MESHFILE)
  +    use arrays,only: dp,elems,elem_cnct,elem_direction,elem_field,&
  +         elem_nodes,elem_ordrs,elem_symmetry,&
  +         nodes,node_xyz,num_elems,&
  +         num_nodes,num_units,units
  +    use indices,only: ne_length,ne_radius,ne_a_A, ne_vol
  +    use other_consts,only: PI
  +    use diagnostics, only: enter_exit
  +    implicit none
  +
  +    character(len=100), intent(in) :: AIRWAY_MESHFILE
  +    ! Local parameters
  +    character(len=100) :: buffer
  +    integer, parameter :: fh = 15
  +    integer :: ios = 0
  +    integer :: line = 0
  +    integer :: i,ibeg,iend,i_ss_end,j,ne,ne0,ne_global,ne_parent,ne_start, &
  +         ngen_parent,np,np0,np_global,&
  +         num_elems_new,num_elems_to_add,num_nodes_new,nunit,nlabel
  +    integer,dimension(1000) :: element_temp,generation, &
  +         parent_element,symmetry_temp
  +    real(dp),dimension(1000) :: length,radius,a_A
  +    character(len=60) :: sub_name
  +
  +    sub_name = 'add_mesh'
  +    call enter_exit(sub_name,1)
  +
  +    open(fh, file=AIRWAY_MESHFILE)
  +
  +    ! ios is negative if an end of record condition is encountered or if
  +    ! an endfile condition was detected.  It is positive if an error was
  +    ! detected.  ios is zero otherwise.
  +    ios=0
  +    line=0
  +    i=0 ! count the number of elements read in
  +
  +    do while (ios == 0)
  +       read(fh, '(A)', iostat=ios) buffer
  +       ! line contains: element, parent element, generation,
  +       !                symmetry, length, outer radius, a/A ratio
  +       ! note that a/A ratio is always 1 for the conducting airways
  +       if (ios == 0) then
  +          line = line + 1
  +          i=i+1
  +          i_ss_end = len(buffer)
  +
  +          do nlabel = 1,7
  +             ibeg = index(buffer," ") + 1 !get location of first integer beyond ws in string
  +             buffer = adjustl(buffer(ibeg:i_ss_end)) ! get info beyond ws, remove leading ws
  +             iend = index(buffer," ") !get location of next ws in sub-string
  +             select case (nlabel)
  +             case (1)
  +                read (buffer(1:iend-1), *, iostat=ios) element_temp(i) ! not used??
  +             case(2)
  +                read (buffer(1:iend-1), *, iostat=ios) parent_element(i)
  +             case(3)
  +                read (buffer(1:iend-1), *, iostat=ios) generation(i)
  +             case(4)
  +                read (buffer(1:iend-1), *, iostat=ios) symmetry_temp(i)
  +             case(5)
  +                read (buffer(1:iend-1), *, iostat=ios) length(i)
  +             case(6)
  +                read (buffer(1:iend-1), *, iostat=ios) radius(i)
  +             case(7)
  +                read (buffer(1:iend-1), *, iostat=ios) a_A(i)
  +             end select
  +          enddo
  +       endif
  +    enddo
  +    close(fh)
  +
  +    num_elems_to_add = i
  +
  +!!! increase the size of node and element arrays to accommodate the additional elements
  +    ! the number of nodes after adding mesh will be:
  +    num_nodes_new = num_nodes + num_units*num_elems_to_add
  +    ! the number of elems after adding mesh will be:
  +    num_elems_new = num_elems + num_units*num_elems_to_add
  +    call reallocate_node_elem_arrays(num_elems_new,num_nodes_new)
  +
  +    ne = num_elems ! the starting local element number
  +    ne_global = elems(ne) ! assumes this is the highest element number (!!!)
  +    np = num_nodes ! the starting local node number
  +    np_global = nodes(np) ! assumes this is the highest node number (!!!)
  +
  +    do nunit = 1,num_units ! for all terminal branches, append the mesh
  +
  +       ne_parent = units(nunit) ! local element number of terminal, to append to
  +       ngen_parent = elem_ordrs(1,ne_parent)
  +       ne_start = ne !starting element number for the unit
  +
  +       do i=1,num_elems_to_add
  +
  +          if(parent_element(i).eq.0)then
  +             ne_parent = units(nunit)
  +          else
  +             ne_parent = ne_start+parent_element(i)
  +          endif
  +
  +          ne0 = ne_parent
  +          np0 = elem_nodes(2,ne0)
  +
  +          ne_global = ne_global + 1 ! new global element number
  +          ne = ne + 1 ! new local element number
  +          np_global = np_global + 1 !new global node number
  +          np = np + 1 ! new local node number
  +
  +          nodes(np) = np_global
  +          elems(ne) = ne_global
  +
  +          elem_nodes(1,ne) = np0
  +          elem_nodes(2,ne) = np
  +
  +          elem_ordrs(1,ne) = ngen_parent + generation(i)
  +          elem_symmetry(ne) = symmetry_temp(i)+1 ! uses 0/1 in file; 1/2 in code
  +
  +          ! record the element connectivity
  +          elem_cnct(-1,0,ne) = 1 ! one parent branch
  +          elem_cnct(-1,1,ne) = ne0 ! store parent element
  +          elem_cnct(1,0,ne0) = elem_cnct(1,0,ne0) + 1
  +          elem_cnct(1,elem_cnct(1,0,ne0),ne0) = ne
  +
  +          ! record the direction and location of the branch
  +          do j=1,3
  +             elem_direction(j,ne) = elem_direction(j,ne0)
  +             node_xyz(j,np) = node_xyz(j,np0) + &
  +                  elem_direction(j,ne)*length(i)
  +          enddo !j
  +
  +          elem_field(ne_length,ne) = length(i)
  +          elem_field(ne_radius,ne) = radius(i)
  +          elem_field(ne_a_A,ne) = a_A(i)
  +          elem_field(ne_vol,ne) = PI*radius(i)**2*length(i)
  +
  +       enddo !i
  +    enddo !nunit
  +
  +    num_nodes = np
  +    num_elems = ne
  +
  +    call element_connectivity_1d
  +    call evaluate_ordering ! calculate new ordering of tree
  +
  +    call enter_exit(sub_name,2)
  +
  +  end subroutine add_mesh
  +!
  +!###################################################################################
  +!
  +!*append_units:* Appends terminal units at the end of a tree structure
  +  subroutine append_units
  +    use arrays,only: elem_cnct,elem_symmetry,elem_units_below,&
  +         num_elems,num_units,units,unit_field
  +    use indices,only: num_nu
  +    use diagnostics, only: enter_exit
  +    implicit none
  +
  +    integer :: ne,ne0,nu
  +    character(len=60) :: sub_name
  +
  +    sub_name = 'append_units'
  +    call enter_exit(sub_name,1)
  +
  +    num_units = 0
  +    DO ne=1,num_elems
  +       IF(elem_cnct(1,0,ne).eq.0)THEN
  +          num_units=num_units+1
  +       ENDIF
  +    ENDDO
  +
  +    if(allocated(units))then !increasing the array size; just overwrite
  +       deallocate(units)
  +       deallocate(unit_field)
  +    endif
  +    allocate(units(num_units))
  +    allocate(unit_field(num_nu,num_units))
  +
  +    unit_field=0.0_dp
  +    units=0
  +    elem_units_below(1:num_elems) = 0 !initialise the number of terminal units below a branch
  +
  +    nu=0
  +    DO ne=1,num_elems
  +       IF(elem_cnct(1,0,ne).eq.0)THEN
  +          nu=nu+1
  +          units(nu)=ne     !Set up units array containing terminals
  +          elem_units_below(ne)=1
  +       ENDIF
  +    ENDDO
  +
  +    ! count the effective number of elements below each branch
  +    do ne=num_elems,2,-1
  +       ne0=elem_cnct(-1,1,ne)
  +       elem_units_below(ne0) = elem_units_below(ne0) &
  +            + elem_units_below(ne)*elem_symmetry(ne)
  +    enddo !ne
  +
  +    call enter_exit(sub_name,2)
  +
  +  end subroutine append_units
  +
  +!
  +!###################################################################################
  +!
  +  subroutine define_1d_elements(ELEMFILE)
  +!*define_1d_elements:* Reads in an element ipelem file to define a geometry
  +    use arrays,only: 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
  +    use diagnostics, only: enter_exit
  +    implicit none
  +
  +    character(LEN=*),intent(in) :: ELEMFILE
  +    !     Local Variables
  +    integer :: ibeg,iend,ierror,i_ss_end,j,ne,ne_global,&
  +         nn,np,np1,np2,np_global
  +    character(LEN=132) :: ctemp1
  +    character(LEN=40) :: sub_string
  +    character(len=60) :: sub_name
  +
  +    sub_name = 'define_1d_elements'
  +    call enter_exit(sub_name,1)
  +
  +    open(10, file=ELEMFILE, status='old')
  +
  +    read_number_of_elements : do
  +       read(unit=10, fmt="(a)", iostat=ierror) ctemp1
  +       if(index(ctemp1, "elements")> 0) then
  +          call get_final_integer(ctemp1,num_elems)
  +          exit read_number_of_elements
  +       endif
  +    end do read_number_of_elements
  +
  +!!! allocate memory for element arrays
  +    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 element arrays
  +    elems=0
  +    elem_nodes=0
  +    elem_symmetry = 1
  +    elem_field = 0.0_dp
  +    expansile = .false.
  +
  +    ne=0
  +
  +    read_an_element : do
  +       !.......read element number
  +       read(unit=10, fmt="(a)", iostat=ierror) ctemp1
  +       if(index(ctemp1, "Element")> 0) then
  +          call get_final_integer(ctemp1,ne_global) !get element number
  +          ne=ne+1
  +          elems(ne)=ne_global
  +
  +          read_element_nodes : do
  +             read(unit=10, fmt="(a)", iostat=ierror) ctemp1
  +             if(index(ctemp1, "global")> 0) then !found the correct line
  +                iend=len(ctemp1)
  +                ibeg=index(ctemp1,":")+1 !get location of first integer in string
  +                sub_string = adjustl(ctemp1(ibeg:iend)) ! get the characters beyond : remove leading blanks
  +                i_ss_end=len(sub_string) !get the end location of the sub-string
  +                ibeg=1
  +                do nn=1,2
  +                   iend=index(sub_string," ") !get location of first blank in sub-string
  +                   read (sub_string(ibeg:iend-1), '(i7)' ) np_global
  +                   call get_local_node(np_global,np) ! get local node np for global node
  +                   elem_nodes(nn,ne)=np ! the local node number, not global
  +                   sub_string = adjustl(sub_string(iend:i_ss_end)) ! get chars beyond blank, remove leading blanks
  +                end do
  +                exit read_element_nodes
  +             endif !index
  +          end do read_element_nodes
  +          if(ne.ge.num_elems) exit read_an_element
  +       endif
  +
  +    end do read_an_element
  +
  +    close(10)
  +
  +    ! 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
  +    enddo
  +
  +    call element_connectivity_1d
  +    call evaluate_ordering
  +
  +    call enter_exit(sub_name,2)
  +
  +  END subroutine define_1d_elements
  +!
  +!###################################################################################
  +!
  +!*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
   !
   !###################################################################################
   !
  
positive.code:     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
  
  !
  !######################################################################################################
  !
  end module geometry
  
  
negative.code:     call enter_exit(sub_name,2)
  
    END subroutine define_rad_from_file
  !
  !##################################################################################
  !
  !*define_rad_from_geom:* Defines vessel or airway radius based on their geometric structure
    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
  