anchor.diff:  !
   !###################################################################################
   !
  -!>volume_of_mesh calculates the volume of an airway mesh including conducting and respiratory airways
  +!*volume_of_mesh:* calculates the volume of an airway mesh including conducting and respiratory airways
     subroutine volume_of_mesh(volume_model,volume_tree)
       use arrays,only: dp,elem_cnct,elem_field,elem_symmetry,&
            num_elems,num_units,units,unit_field
  
positive.code:     call enter_exit(sub_name,2)
  
    END subroutine define_rad_from_geom
  
  
  
  !inputs element list and radius scheme
  
  
  !
  !###################################################################################
  !
  !>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
  
  !
  !###########################################################################
  !
  !>element_connectivity_1d  Calculates element connectivity in 1D and stores in elelem_cnct
    subroutine element_connectivity_1d
      use arrays,only: elem_cnct,elem_nodes,elems_at_node,num_elems,num_nodes
      use diagnostics, only: enter_exit
  
negative.code: 
    !Interfaces
    private
    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
  !
  !###################################################################################
  !
  