anchor.diff:    if ( ios /= 0 ) stop "Error opening file 21"
   
     do ii = 1, num_nodes
  -    write(21,*) xcoordsN(ii), GlobalX(ii)
  +    write(21,*) xcoords(ii), GlobalX(ii)
     end do
     close(unit=21, iostat=ios)
     if ( ios /= 0 ) stop "Error closing file unit 21"
   
  -  if ( allocated(orderN) )    deallocate(orderN)
  -  if ( allocated(elem_connN) ) deallocate(elem_connN)
  -  if ( allocated(xcoordsN) )  deallocate(xcoordsN)
  +  if ( allocated(order) )    deallocate(order)
  +  if ( allocated(elem_conn) ) deallocate(elem_conn)
  +  if ( allocated(xcoords) )  deallocate(xcoords)
     if ( allocated(GlobalA) )  deallocate(GlobalA)
     if ( allocated(GlobalB) )  deallocate(GlobalB)
     if ( allocated(GlobalX) )  deallocate(GlobalX)
  
positive.code:   integer             :: ii, num_nodes, ios
    real(wp), parameter :: diff = 0.1_wp, vel = -1.0_wp
  
    integer,  dimension(:),   allocatable :: orderN
    integer,  dimension(:,:), allocatable :: elem_connN
    real(wp), dimension(:),   allocatable :: xcoordsN, GlobalB, GlobalX
    real(wp), dimension(:,:), allocatable :: IeN, GlobalA
    character(len=50) :: arg
  
    call get_command_argument_wrapper(arg)
    call read_gmsh_file_1D(arg, num_nodes, orderN, elem_connN, xcoordsN)
  
    allocate(GlobalA(num_nodes, num_nodes))
    allocate(GlobalB(num_nodes))
  
negative.code:   call r8mat_print(num_nodes, 1, GlobalB, 'Global RHS:')
  
    ! Add elemental stiffness matrices to Global Stiffness Matrix
    do ii = 1, size(order)
      ! Reallocate elemental stiffness matrix
      allocate(Ie(order(ii)+1, order(ii)+1))
      call getIe(1, 1, xcoords(elem_conn(ii,:)), Ie)
      ! call r8mat_print(order(ii)+1, order(ii)+1, Ie, 'Elemental Stiffness Matrix:')
  
      GlobalA(elem_conn(ii,:), elem_conn(ii,:)) = &
          GlobalA(elem_conn(ii,:), elem_conn(ii,:)) - diff*Ie
  
      call getIe(0, 1, xcoords(elem_conn(ii,:)), Ie)
      ! call r8mat_print(order(ii)+1, order(ii)+1, Ie, 'Elemental Stiffness Matrix:')
  
      GlobalA(elem_conn(ii,:), elem_conn(ii,:)) = &
          GlobalA(elem_conn(ii,:), elem_conn(ii,:)) - vel*Ie
  
      ! Deallocate elemental stiffness matrix after every loop
      deallocate(Ie)
      ! stop
    end do
  
  