MODULE Lattice_Geometry
      USE Precision
      USE Error_Handling
      USE System_Monitors
      USE Random_Numbers
      USE Graph_Algorithms
      USE Network_Data_Structures
      USE Network_Graphics
      
      IMPLICIT NONE
      PUBLIC::InitializeLattice,CreateLatticeNetwork,DestroyLatticeNetwork
      PRIVATE
      
      INTEGER,PARAMETER,PUBLIC::n_directions=2
      INTEGER,DIMENSION(n_dim),PUBLIC,SAVE::lengths
      CHARACTER,PUBLIC,SAVE::component_extraction='C'
      INTEGER,DIMENSION(n_dim),PUBLIC,SAVE::&
      source_layer_thickness=1,sink_layer_thickness=1,&
      injection_layer_thickness=1,outflow_layer_thickness=1
      INTEGER,DIMENSION(n_dim,n_directions),PARAMETER::displacements=&
      RESHAPE(SOURCE=(/1,0,0,1/),SHAPE=(/n_dim,n_directions/))
      INTEGER,PARAMETER,PUBLIC::bc_periodic=1,bc_injection=2,bc_plate=3,bc_source_sink=4
      LOGICAL,DIMENSION(4,n_dim),PUBLIC,SAVE::boundaries=.FALSE.
      INTEGER(KIND=i_wp),DIMENSION(n_dim),PUBLIC::sources,sinks
      INTEGER(KIND=i_wp),PUBLIC,SAVE::n_injection_nodes,n_outflow_nodes
      REAL,PUBLIC,SAVE::node_dilution=-1.0,arc_dilution=-1.0
      INTEGER(KIND=i_byte),DIMENSION(:),ALLOCATABLE,PUBLIC::arcs_status,nodes_status
      INTEGER(KIND=i_byte),PARAMETER,PUBLIC::empty_node=0,interior_node=1,&
      plate_node=2,to_source_node=3,to_sink_node=4,&
      injection_node=5,outflow_node=6,source_node=7,sink_node=8
      INTEGER(KIND=i_byte),PARAMETER,PUBLIC::empty_arc=0,interior_arc=1,periodic_bc_arc=2,&
      plate_arc=3,source_arc=4,sink_arc=5
      
      NAMELIST/LatticeOptions/lengths,component_extraction,source_layer_thickness,sink_layer_thickness,injection_layer_thickness,o&
     &utflow_layer_thickness,node_dilution,arc_dilution,boundaries
      
      CONTAINS
      
      
      SUBROUTINE InitializeLattice()
      IMPLICIT NONE
      
      REWIND(UNIT=program_options_unit)
      READ(UNIT=program_options_unit,NML=LatticeOptions,IOSTAT=error_status)
      IF(error_status/=0)&
      CALL CriticalError(message="NAMELIST LatticeOptions was not read successfully"//&
      " from file "//TRIM(options_file),caller="InitializeLattice")
      WRITE(UNIT=message_log_unit,NML=LatticeOptions)
      
      END SUBROUTINE InitializeLattice
      
      
      
      SUBROUTINE CreateLatticeNetwork()
      IMPLICIT NONE
      
      INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::points_count,nodes_labels,nodes_count
      
      
      INTEGER(KIND=i_wp),DIMENSION(:,:),ALLOCATABLE::heads_tails_buffer
      REAL,DIMENSION(:,:),ALLOCATABLE::coords_buffer
      INTEGER(KIND=i_byte),DIMENSION(:),ALLOCATABLE::status_buffer
      
      INTEGER(KIND=i_wp)::n_points,n_links,n_special_points,n_special_links,&
      n_selected_nodes,n_selected_arcs,n_selected_special_arcs
      INTEGER(KIND=i_wp)::link,point,node,arc,neighbor_point,node_counter,arc_counter,&
      point_counter,point_temp,selected_label,selected_size
      INTEGER::dim,direction,dim_temp
      INTEGER(KIND=i_byte)::arc_status,node_status
      INTEGER(KIND=i_wp),DIMENSION(n_dim)::lengths_product
      INTEGER,DIMENSION(n_dim)::neighbor_coords,point_coords
      REAL::dice,real_temp
      LOGICAL::dilute_nodes,dilute_arcs,extraction_success,&
      node_on_plate,special_arc,renumber_nodes
      
      
      lengths_product(1)=1
      DO dim=2,n_dim
      lengths_product(dim)=lengths_product(dim-1)*lengths(dim-1)
      END DO
      n_points=PRODUCT(lengths)
      n_links=n_points*n_directions
      sources=0;
      sinks=0;
      
      
      n_special_points=0;
      n_special_links=0;
      
      DO dim=1,n_dim
      PlateArcs:IF(boundaries(bc_plate,dim))THEN
      DO direction=1,n_directions
      IF(displacements(dim,direction)==0)THEN
      n_special_links=n_special_links+2*n_points/lengths(dim)
      n_links=n_links-2*n_points/lengths(dim)
      END IF
      END DO
      END IF PlateArcs
      SourcesSinks:IF(boundaries(bc_source_sink,dim))THEN
      IF(source_layer_thickness(dim)>0)THEN
      n_special_links=n_special_links+&
      source_layer_thickness(dim)*n_points/lengths(dim)
      n_special_points=n_special_points+1
      END IF
      IF(sink_layer_thickness(dim)>0)THEN
      n_special_links=n_special_links+&
      sink_layer_thickness(dim)*n_points/lengths(dim)
      n_special_points=n_special_points+1
      END IF
      END IF SourcesSinks
      END DO
      
      
      ALLOCATE(nodes_coords(n_dim,-n_special_points:n_points),STAT=error_status)
      CALL RecordAllocation(n_elements=n_dim*(n_points+n_special_points+1),mold=1.0,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      nodes_coords(1:n_dim,0)=0.0
      ALLOCATE(nodes_status(-n_special_points:n_points),STAT=error_status)
      CALL RecordAllocation(n_elements=(n_points+n_special_points+1),mold=1_i_byte,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      ALLOCATE(heads_tails(2,-n_special_links:n_links),STAT=error_status)
      CALL RecordAllocation(n_elements=2*(n_links+n_special_links+1),mold=1_i_wp,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      heads_tails(1:2,0)=0
      ALLOCATE(arcs_status(-n_special_links:n_links),STAT=error_status)
      CALL RecordAllocation(n_elements=(n_links+n_special_links+1),mold=1_i_byte,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      nodes_status=empty_node;
      arcs_status=empty_arc
      
      
      
      dilute_nodes=((node_dilution<1.0).AND.(node_dilution>0.0))
      dilute_arcs=((arc_dilution<1.0).AND.(arc_dilution>0.0))
      
      
      IF(dilute_nodes)THEN
      ALLOCATE(points_count(n_points),STAT=error_status)
      CALL RecordAllocation(n_elements=n_points,mold=1_i_wp,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      points_count=0
      point_counter=0
      DiluteNodes:DO point=1,n_points
      node_on_plate=.FALSE.
      
      point_temp=point-1
      DO dim_temp=n_dim,1,-1
      point_coords(dim_temp)=point_temp/lengths_product(dim_temp)+1
      point_temp=point_temp-(point_coords(dim_temp)-1)*lengths_product(dim_temp)
      END DO
      PlateNodes:DO dim=1,n_dim
      IF(boundaries(bc_plate,dim))THEN
      IF(point_coords(dim)==1.OR.point_coords(dim)==lengths(dim))&
      node_on_plate=.TRUE.
      END IF
      END DO PlateNodes
      IF(.NOT.node_on_plate)THEN
      CALL RandomUniform(dice)
      IF(dice>node_dilution)CYCLE DiluteNodes
      END IF
      point_counter=point_counter+1
      points_count(point)=point_counter
      END DO DiluteNodes
      END IF
      
      n_nodes=0;
      n_arcs=0;
      n_special_nodes=0;
      n_special_arcs=0;
      
      
      SourceSinkNodes:DO dim=1,n_dim
      IF(boundaries(bc_source_sink,dim))THEN
      IF(source_layer_thickness(dim)>0)THEN
      n_special_nodes=n_special_nodes+1
      sources(dim)=-n_special_nodes
      nodes_status(-n_special_nodes)=source_node
      nodes_coords(1:n_dim,-n_special_nodes)=0.5*REAL(lengths+1)
      nodes_coords(dim,-n_special_nodes)=1.0-0.1*REAL(lengths(dim))
      END IF
      IF(sink_layer_thickness(dim)>0)THEN
      n_special_nodes=n_special_nodes+1
      sinks(dim)=-n_special_nodes
      nodes_status(-n_special_nodes)=sink_node
      nodes_coords(1:n_dim,-n_special_nodes)=0.5*REAL(lengths+1)
      nodes_coords(dim,-n_special_nodes)=1.1*REAL(lengths(dim))
      END IF
      END IF
      END DO SourceSinkNodes
      
      
      ForAllPoints:DO point=1,n_points
      node_status=interior_node
      
      IF(dilute_nodes)THEN
      IF(points_count(point)==0)CYCLE ForAllPoints
      n_nodes=points_count(point)
      ELSE
      n_nodes=point
      END IF
      
      
      point_temp=point-1
      DO dim_temp=n_dim,1,-1
      point_coords(dim_temp)=point_temp/lengths_product(dim_temp)+1
      point_temp=point_temp-(point_coords(dim_temp)-1)*lengths_product(dim_temp)
      END DO
      nodes_coords(1:n_dim,n_nodes)=REAL(point_coords)
      
      PlateOrInjectionNodes:DO dim=1,n_dim
      IF(boundaries(bc_plate,dim))THEN
      IF(point_coords(dim)==1.OR.point_coords(dim)==lengths(dim))THEN
      node_status=MAX(node_status,plate_node)
      END IF
      END IF
      IF(boundaries(bc_injection,dim))THEN
      IF(point_coords(dim)<=injection_layer_thickness(dim))THEN
      node_status=MAX(node_status,injection_node)
      ELSE IF(point_coords(dim)>lengths(dim)-outflow_layer_thickness(dim))THEN
      node_status=MAX(node_status,outflow_node)
      END IF
      END IF
      IF(boundaries(bc_source_sink,dim))THEN
      IF(point_coords(dim)<=source_layer_thickness(dim))THEN
      node_status=MAX(node_status,to_source_node)
      n_special_arcs=n_special_arcs+1
      arcs_status(-n_special_arcs)=MAX(arcs_status(-n_special_arcs),source_arc)
      IF(dilute_nodes)THEN
      heads_tails(1:2,-n_special_arcs)=(/sources(dim),points_count(point)/)
      ELSE
      heads_tails(1:2,-n_special_arcs)=(/sources(dim),point/)
      END IF
      ELSE IF(point_coords(dim)>lengths(dim)-sink_layer_thickness(dim))THEN
      node_status=MAX(node_status,to_sink_node)
      n_special_arcs=n_special_arcs+1
      arcs_status(-n_special_arcs)=MAX(arcs_status(-n_special_arcs),sink_arc)
      IF(dilute_nodes)THEN
      heads_tails(1:2,-n_special_arcs)=(/points_count(point),sinks(dim)/)
      ELSE
      heads_tails(1:2,-n_special_arcs)=(/point,sinks(dim)/)
      END IF
      END IF
      END IF
      END DO PlateOrInjectionNodes
      nodes_status(n_nodes)=node_status
      
      
      ForAllDirections:DO direction=1,n_directions
      arc_status=interior_arc
      special_arc=.FALSE.
      
      neighbor_coords=point_coords+displacements(:,direction)
      
      Plates:DO dim=1,n_dim
      BoundaryPlate:IF(boundaries(bc_plate,dim))THEN
      IF(((point_coords(dim)==lengths(dim)).AND.(neighbor_coords(dim)==lengths(dim)))&
      .OR.((point_coords(dim)==1).AND.(neighbor_coords(dim)==1)))THEN
      special_arc=.TRUE.
      arc_status=MAX(plate_arc,arc_status)
      END IF
      END IF BoundaryPlate
      END DO Plates
      
      IF(.NOT.special_arc.AND.dilute_arcs)THEN
      CALL RandomUniform(dice)
      IF(dice>arc_dilution)CYCLE ForAllDirections
      END IF
      
      Periodics:DO dim=1,n_dim
      BoundaryCrossing:IF((neighbor_coords(dim)>lengths(dim))&
      .OR.(neighbor_coords(dim)<1))THEN
      IF(boundaries(bc_periodic,dim))THEN
      IF(neighbor_coords(dim)>lengths(dim))&
      neighbor_coords(dim)=MOD(neighbor_coords(dim),lengths(dim))
      IF(neighbor_coords(dim)<1)&
      neighbor_coords(dim)=lengths(dim)+MOD(neighbor_coords(dim),lengths(dim))
      arc_status=MAX(periodic_bc_arc,arc_status)
      ELSE
      CYCLE ForAllDirections
      END IF
      END IF BoundaryCrossing
      END DO Periodics
      
      
      neighbor_point=1
      DO dim_temp=1,n_dim
      neighbor_point=neighbor_point+(neighbor_coords(dim_temp)-1)*lengths_product(dim_temp)
      END DO
      
      IF(dilute_nodes)THEN
      IF(points_count(neighbor_point)==0)CYCLE ForAllDirections
      IF(.NOT.special_arc)THEN
      n_arcs=n_arcs+1
      heads_tails(1:2,n_arcs)=(/points_count(point),points_count(neighbor_point)/)
      arcs_status(n_arcs)=arc_status
      ELSE
      n_special_arcs=n_special_arcs+1
      heads_tails(1:2,-n_special_arcs)=(/points_count(point),points_count(neighbor_point)/)
      arcs_status(-n_special_arcs)=arc_status
      END IF
      ELSE
      IF(.NOT.special_arc)THEN
      n_arcs=n_arcs+1
      heads_tails(1:2,n_arcs)=(/point,neighbor_point/)
      arcs_status(n_arcs)=arc_status
      ELSE
      n_special_arcs=n_special_arcs+1
      heads_tails(1:2,-n_special_arcs)=(/point,neighbor_point/)
      arcs_status(-n_special_arcs)=arc_status
      END IF
      END IF
      
      END DO ForAllDirections
      
      END DO ForAllPoints
      IF(dilute_nodes)THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(points_count),KIND=i_wp),&
      mold=1_i_wp,caller="CreateLatticeNetwork")
      DEALLOCATE(points_count)
      END IF
      
      
      
      n_injection_nodes=0
      n_outflow_nodes=0
      
      ExtractSelectedComponent:IF(dilute_nodes.OR.dilute_arcs)THEN
      
      WRITE(message_log_unit,*)"Starting connected-component extraction with:"
      WRITE(message_log_unit,*)"n_arcs=",n_arcs," n_nodes=",n_nodes,&
      " n_special_nodes=",n_special_nodes," n_special_arcs=",n_special_arcs
      
      ALLOCATE(nodes_mask(-n_special_nodes:n_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_nodes+n_nodes+1,mold=.TRUE._l_wp,&
      caller="CreateNetwork",alloc_status=error_status)
      ALLOCATE(arcs_mask(-n_special_arcs:n_arcs))
      CALL RecordAllocation(n_elements=n_special_arcs+n_arcs+1,mold=.TRUE._l_wp,&
      caller="CreateNetwork",alloc_status=error_status)
      
      ALLOCATE(nodes_labels(-n_special_nodes:n_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_nodes+n_nodes+1,mold=1_i_wp,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      nodes_labels(0)=0
      
      SELECT CASE(component_extraction)
      CASE('b','B')
      CALL BackboneCycles(heads_tails=heads_tails(:,-n_special_arcs:n_arcs),&
      node_offset=n_special_nodes,supernodes=nodes_labels,&
      largest_supernode=selected_label,&
      largest_supernode_size=selected_size)
      WRITE(message_log_unit,*)"The largest backbone component has: ",&
      selected_size," nodes"
      CASE DEFAULT
      CALL ConnectedComponents(heads_tails=heads_tails(:,-n_special_arcs:n_arcs),&
      node_offset=n_special_nodes,labels=nodes_labels,&
      largest_component_label=selected_label,&
      largest_component_size=selected_size)
      WRITE(message_log_unit,*)"The largest connected component has: ",&
      selected_size," nodes"
      ENDSELECT
      
      n_selected_nodes=0
      n_selected_arcs=0
      n_selected_special_arcs=0
      extraction_success=.TRUE.
      
      
      MaskSelectedNodes:DO node=-n_special_nodes,n_nodes
      SelectNodes:IF(nodes_labels(node)==selected_label)THEN
      nodes_mask(node)=.TRUE._l_wp
      IF(node>0)n_selected_nodes=n_selected_nodes+1
      IF(nodes_status(node)==injection_node)THEN
      n_injection_nodes=n_injection_nodes+1
      ELSE IF(nodes_status(node)==outflow_node)THEN
      n_outflow_nodes=n_outflow_nodes+1
      END IF
      ELSE
      IF((node<0).OR.(nodes_status(node)==plate_node))THEN
      extraction_success=.FALSE.
      EXIT MaskSelectedNodes
      END IF
      nodes_mask(node)=.FALSE._l_wp
      END IF SelectNodes
      END DO MaskSelectedNodes
      
      CheckConnectedness:IF(.NOT.extraction_success)THEN
      CALL CriticalError(message="The sources/sinks and/or plate nodes were not in the largest connected component",caller="Create&
     &LatticeNetwork")
      ELSE
      DO dim=1,n_dim
      IF(boundaries(bc_injection,dim))THEN
      IF((injection_layer_thickness(dim)>0).AND.(n_injection_nodes==0))&
      CALL CriticalError(message="No injection sites in the largest component",caller="CreateLatticeNetwork")
      IF((outflow_layer_thickness(dim)>0).AND.(n_outflow_nodes==0))&
      CALL CriticalError(message="No outflow sites in the largest component",caller="CreateLatticeNetwork")
      END IF
      END DO
      END IF CheckConnectedness
      
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_labels),KIND=i_wp),&
      mold=1_i_wp,caller="CreateLatticeNetwork")
      DEALLOCATE(nodes_labels)
      
      MaskSelectedSpecialArcs:DO arc=-n_special_arcs,-1
      IF(nodes_mask(heads_tails(1,arc)).AND.nodes_mask(heads_tails(2,arc)))THEN
      n_selected_special_arcs=n_selected_special_arcs+1
      arcs_mask(arc)=.TRUE._l_wp
      ELSE
      arcs_mask(arc)=.FALSE._l_wp
      END IF
      END DO MaskSelectedSpecialArcs
      
      MaskSelectedArcs:DO arc=1,n_arcs
      IF(nodes_mask(heads_tails(1,arc)).AND.nodes_mask(heads_tails(2,arc)))THEN
      n_selected_arcs=n_selected_arcs+1
      arcs_mask(arc)=.TRUE._l_wp
      ELSE
      arcs_mask(arc)=.FALSE._l_wp
      END IF
      END DO MaskSelectedArcs
      
      
      renumber_nodes=.FALSE.
      IF((n_selected_nodes/=n_points).OR.(n_special_nodes/=n_special_points))THEN
      WRITE(message_log_unit,*)"Reallocating node arrays..."
      
      ALLOCATE(coords_buffer(n_dim,-n_special_nodes:n_selected_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_dim*(n_special_nodes+n_selected_nodes+1),mold=1.0,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      ALLOCATE(status_buffer(-n_special_nodes:n_selected_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_nodes+n_selected_nodes+1,mold=1_i_byte,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      coords_buffer(1:n_dim,0)=0.0
      status_buffer(0)=empty_node
      
      
      IF(n_selected_nodes/=n_nodes)THEN
      renumber_nodes=.TRUE.
      WRITE(message_log_unit,*)"Repacking node arrays..."
      ALLOCATE(nodes_count(-n_special_nodes:n_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_nodes+n_nodes+1,mold=1_i_wp,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      
      nodes_count(0)=0
      DO node=-n_special_nodes,-1
      nodes_count(node)=node
      END DO
      node_counter=0
      DO node=1,n_nodes
      IF(nodes_mask(node))THEN
      node_counter=node_counter+1
      nodes_count(node)=node_counter
      END IF
      END DO
      
      
      
      coords_buffer(1:n_dim,:-1)=nodes_coords(1:n_dim,:-1)
      status_buffer(:-1)=nodes_status(:-1)
      node_counter=0
      PackNodes:DO node=1,n_nodes
      IF(nodes_mask(node))THEN
      node_counter=node_counter+1
      coords_buffer(1:n_dim,node_counter)=nodes_coords(1:n_dim,node)
      status_buffer(node_counter)=nodes_status(node)
      END IF
      END DO PackNodes
      
      
      ELSE
      
      coords_buffer=nodes_coords(1:n_dim,-n_special_nodes:n_selected_nodes)
      status_buffer=nodes_status(-n_special_nodes:n_selected_nodes)
      
      
      END IF
      
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_coords),KIND=i_wp),&
      mold=1.0,caller="CreateLatticeNetwork")
      DEALLOCATE(nodes_coords)
      ALLOCATE(nodes_coords(n_dim,-n_special_nodes:n_selected_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_dim*(n_special_nodes+n_selected_nodes+1),mold=1.0,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      nodes_coords=coords_buffer
      CALL RecordAllocation(n_elements=-INT(SIZE(coords_buffer),KIND=i_wp),&
      mold=1.0,caller="CreateLatticeNetwork")
      DEALLOCATE(coords_buffer)
      
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_status),KIND=i_wp),&
      mold=1_i_byte,caller="CreateLatticeNetwork")
      DEALLOCATE(nodes_status)
      ALLOCATE(nodes_status(-n_special_nodes:n_selected_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_nodes+n_selected_nodes+1,mold=1_i_byte,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      nodes_status=status_buffer
      CALL RecordAllocation(n_elements=-INT(SIZE(status_buffer),KIND=i_wp),&
      mold=1_i_byte,caller="CreateLatticeNetwork")
      DEALLOCATE(status_buffer)
      
      
      END IF
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_mask),KIND=i_wp),&
      mold=.TRUE._l_wp,caller="CreateLatticeNetwork")
      DEALLOCATE(nodes_mask)
      
      IF((n_selected_arcs/=n_links).OR.(n_selected_special_arcs/=n_special_links))THEN
      
      WRITE(message_log_unit,*)"Reallocating arc arrays..."
      
      
      ALLOCATE(heads_tails_buffer(2,-n_selected_special_arcs:n_selected_arcs),&
      STAT=error_status)
      CALL RecordAllocation(n_elements=2*(n_selected_arcs+n_selected_special_arcs+1),&
      mold=1_i_wp,caller="CreateLatticeNetwork",alloc_status=error_status)
      ALLOCATE(status_buffer(-n_selected_special_arcs:n_selected_arcs),STAT=error_status)
      CALL RecordAllocation(n_elements=n_selected_arcs+n_selected_special_arcs+1,mold=1_i_byte,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      heads_tails_buffer(1:2,0)=0
      status_buffer(0)=empty_arc
      
      
      
      IF((n_selected_arcs/=n_arcs).OR.(n_selected_special_arcs/=n_special_arcs))THEN
      
      WRITE(message_log_unit,*)"Repacking arc arrays..."
      IF(renumber_nodes)THEN
      
      arc_counter=0
      DO arc=-1,-n_special_arcs,-1
      IF(arcs_mask(arc))THEN
      arc_counter=arc_counter-1
      heads_tails_buffer(1:2,arc_counter)=nodes_count(heads_tails(1:2,arc))
      status_buffer(arc_counter)=arcs_status(arc)
      END IF
      END DO
      
      arc_counter=0
      PackArcsRenumbered:DO arc=1,n_arcs
      IF(arcs_mask(arc))THEN
      arc_counter=arc_counter+1
      heads_tails_buffer(1:2,arc_counter)=nodes_count(heads_tails(1:2,arc))
      status_buffer(arc_counter)=arcs_status(arc)
      END IF
      END DO PackArcsRenumbered
      
      
      ELSE
      
      arc_counter=0
      DO arc=-1,-n_special_arcs,-1
      IF(arcs_mask(arc))THEN
      arc_counter=arc_counter-1
      heads_tails_buffer(1:2,arc_counter)=heads_tails(1:2,arc)
      status_buffer(arc_counter)=arcs_status(arc)
      END IF
      END DO
      
      arc_counter=0
      PackArcs:DO arc=1,n_arcs
      IF(arcs_mask(arc))THEN
      arc_counter=arc_counter+1
      heads_tails_buffer(1:2,arc_counter)=heads_tails(1:2,arc)
      status_buffer(arc_counter)=arcs_status(arc)
      END IF
      END DO PackArcs
      
      
      END IF
      ELSE
      IF(renumber_nodes)THEN
      
      DO arc=-n_selected_special_arcs,n_selected_arcs
      heads_tails_buffer(1:2,arc)=nodes_count(heads_tails(1:2,arc))
      END DO
      status_buffer=arcs_status(-n_selected_special_arcs:n_selected_arcs)
      
      
      ELSE
      
      heads_tails_buffer=heads_tails(1:2,-n_selected_special_arcs:n_selected_arcs)
      status_buffer=arcs_status(-n_selected_special_arcs:n_selected_arcs)
      
      
      END IF
      END IF
      
      CALL RecordAllocation(n_elements=-INT(SIZE(heads_tails),KIND=i_wp),&
      mold=1_i_wp,caller="CreateLatticeNetwork")
      DEALLOCATE(heads_tails)
      ALLOCATE(heads_tails(2,-n_selected_special_arcs:n_selected_arcs),&
      STAT=error_status)
      CALL RecordAllocation(n_elements=2*(n_selected_arcs+n_selected_special_arcs+1),&
      mold=1_i_wp,caller="CreateLatticeNetwork",alloc_status=error_status)
      heads_tails=heads_tails_buffer
      CALL RecordAllocation(n_elements=-INT(SIZE(heads_tails_buffer),KIND=i_wp),&
      mold=1_i_wp,caller="CreateLatticeNetwork")
      DEALLOCATE(heads_tails_buffer)
      
      CALL RecordAllocation(n_elements=-INT(SIZE(arcs_status),KIND=i_wp),&
      mold=1_i_byte,caller="CreateLatticeNetwork")
      DEALLOCATE(arcs_status)
      ALLOCATE(arcs_status(-n_selected_special_arcs:n_selected_arcs),&
      STAT=error_status)
      CALL RecordAllocation(n_elements=n_selected_arcs+n_selected_special_arcs+1,mold=1_i_byte,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      arcs_status=status_buffer
      CALL RecordAllocation(n_elements=-INT(SIZE(status_buffer),KIND=i_wp),&
      mold=1_i_byte,caller="CreateLatticeNetwork")
      DEALLOCATE(status_buffer)
      
      
      END IF
      IF(renumber_nodes)THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_count),KIND=i_wp),&
      mold=1_i_wp,caller="CreateLatticeNetwork")
      DEALLOCATE(nodes_count)
      END IF
      CALL RecordAllocation(n_elements=-INT(SIZE(arcs_mask),KIND=i_wp),&
      mold=.TRUE._l_wp,caller="CreateLatticeNetwork")
      DEALLOCATE(arcs_mask)
      
      n_nodes=n_selected_nodes
      n_arcs=n_selected_arcs
      n_special_arcs=n_selected_special_arcs
      
      ELSE
      
      n_selected_nodes=n_nodes
      n_selected_arcs=n_arcs
      n_selected_special_arcs=n_special_arcs
      
      IF((n_nodes/=n_points).OR.(n_special_nodes/=n_special_points))THEN
      WRITE(message_log_unit,*)"Reallocating node arrays..."
      
      ALLOCATE(coords_buffer(n_dim,-n_special_nodes:n_selected_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_dim*(n_special_nodes+n_selected_nodes+1),mold=1.0,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      ALLOCATE(status_buffer(-n_special_nodes:n_selected_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_nodes+n_selected_nodes+1,mold=1_i_byte,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      coords_buffer(1:n_dim,0)=0.0
      status_buffer(0)=empty_node
      
      
      
      coords_buffer=nodes_coords(1:n_dim,-n_special_nodes:n_selected_nodes)
      status_buffer=nodes_status(-n_special_nodes:n_selected_nodes)
      
      
      
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_coords),KIND=i_wp),&
      mold=1.0,caller="CreateLatticeNetwork")
      DEALLOCATE(nodes_coords)
      ALLOCATE(nodes_coords(n_dim,-n_special_nodes:n_selected_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_dim*(n_special_nodes+n_selected_nodes+1),mold=1.0,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      nodes_coords=coords_buffer
      CALL RecordAllocation(n_elements=-INT(SIZE(coords_buffer),KIND=i_wp),&
      mold=1.0,caller="CreateLatticeNetwork")
      DEALLOCATE(coords_buffer)
      
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_status),KIND=i_wp),&
      mold=1_i_byte,caller="CreateLatticeNetwork")
      DEALLOCATE(nodes_status)
      ALLOCATE(nodes_status(-n_special_nodes:n_selected_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_nodes+n_selected_nodes+1,mold=1_i_byte,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      nodes_status=status_buffer
      CALL RecordAllocation(n_elements=-INT(SIZE(status_buffer),KIND=i_wp),&
      mold=1_i_byte,caller="CreateLatticeNetwork")
      DEALLOCATE(status_buffer)
      
      
      END IF
      
      IF((n_arcs/=n_links).OR.(n_special_arcs/=n_special_links))THEN
      WRITE(message_log_unit,*)"Reallocating arc arrays..."
      
      ALLOCATE(heads_tails_buffer(2,-n_selected_special_arcs:n_selected_arcs),&
      STAT=error_status)
      CALL RecordAllocation(n_elements=2*(n_selected_arcs+n_selected_special_arcs+1),&
      mold=1_i_wp,caller="CreateLatticeNetwork",alloc_status=error_status)
      ALLOCATE(status_buffer(-n_selected_special_arcs:n_selected_arcs),STAT=error_status)
      CALL RecordAllocation(n_elements=n_selected_arcs+n_selected_special_arcs+1,mold=1_i_byte,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      heads_tails_buffer(1:2,0)=0
      status_buffer(0)=empty_arc
      
      
      
      heads_tails_buffer=heads_tails(1:2,-n_selected_special_arcs:n_selected_arcs)
      status_buffer=arcs_status(-n_selected_special_arcs:n_selected_arcs)
      
      
      
      CALL RecordAllocation(n_elements=-INT(SIZE(heads_tails),KIND=i_wp),&
      mold=1_i_wp,caller="CreateLatticeNetwork")
      DEALLOCATE(heads_tails)
      ALLOCATE(heads_tails(2,-n_selected_special_arcs:n_selected_arcs),&
      STAT=error_status)
      CALL RecordAllocation(n_elements=2*(n_selected_arcs+n_selected_special_arcs+1),&
      mold=1_i_wp,caller="CreateLatticeNetwork",alloc_status=error_status)
      heads_tails=heads_tails_buffer
      CALL RecordAllocation(n_elements=-INT(SIZE(heads_tails_buffer),KIND=i_wp),&
      mold=1_i_wp,caller="CreateLatticeNetwork")
      DEALLOCATE(heads_tails_buffer)
      
      CALL RecordAllocation(n_elements=-INT(SIZE(arcs_status),KIND=i_wp),&
      mold=1_i_byte,caller="CreateLatticeNetwork")
      DEALLOCATE(arcs_status)
      ALLOCATE(arcs_status(-n_selected_special_arcs:n_selected_arcs),&
      STAT=error_status)
      CALL RecordAllocation(n_elements=n_selected_arcs+n_selected_special_arcs+1,mold=1_i_byte,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      arcs_status=status_buffer
      CALL RecordAllocation(n_elements=-INT(SIZE(status_buffer),KIND=i_wp),&
      mold=1_i_byte,caller="CreateLatticeNetwork")
      DEALLOCATE(status_buffer)
      
      
      END IF
      
      IF(ANY(boundaries(bc_injection,:)))THEN
      DO node=1,n_nodes
      IF(nodes_status(node)==injection_node)THEN
      n_injection_nodes=n_injection_nodes+1
      ELSE IF(nodes_status(node)==outflow_node)THEN
      n_outflow_nodes=n_outflow_nodes+1
      END IF
      END DO
      END IF
      
      END IF ExtractSelectedComponent
      
      IF((n_injection_nodes/=0).OR.(n_outflow_nodes/=0))&
      WRITE(message_log_unit,*)"There are ",n_injection_nodes," injection sites, and ",&
      n_outflow_nodes," outflow sites in the network."
      
      
      
      END SUBROUTINE CreateLatticeNetwork
      
      
      
      SUBROUTINE DestroyLatticeNetwork()
      IMPLICIT NONE
      
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_coords),KIND=i_wp),&
      mold=1.0,caller="DestroyLatticeNetwork")
      DEALLOCATE(nodes_coords)
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_status),KIND=i_wp),&
      mold=1_i_byte,caller="DestroyLatticeNetwork")
      DEALLOCATE(nodes_status)
      CALL RecordAllocation(n_elements=-INT(SIZE(heads_tails),KIND=i_wp),&
      mold=1_i_wp,caller="DestroyLatticeNetwork")
      DEALLOCATE(heads_tails)
      CALL RecordAllocation(n_elements=-INT(SIZE(arcs_status),KIND=i_wp),&
      mold=1_i_byte,caller="DestroyLatticeNetwork")
      DEALLOCATE(arcs_status)
      
      END SUBROUTINE DestroyLatticeNetwork
      
      
      
      END MODULE Lattice_Geometry