How to combine the largest subtours together

Here is a somewhat modified version of the subtour elimination subroutine that works with the largest subtour for defining the constraints:

procedure remove_subtours
  !declarations of the procedure
  declarations
   ROUTE,LARGEST,ALLRETAIL: set of integer
  end-declarations
 
!Determining the next retail shop visited in current solution
  forall(i in RETAIL) 
   NEXTRETAIL(i):= integer(round(getsol(sum(j in RETAIL) j*ROAD(i,j) )))
 
! Print the solution of each tour with the other procedure
  print_allsubsolution
  
  
! Subtour which contains the starting point 0  
  ROUTE:={}
  start:=0
  repeat
   ROUTE+={start}
   start:= NEXTRETAIL(start)
  until start=0
  size:=getsize(ROUTE)
 
! Find the biggest subtour
  if size < NRETAIL then
    LARGEST:=ROUTE
    ALLRETAIL:=ROUTE
    forall(i in RETAIL) do
     if(i not in ALLRETAIL) then
      ROUTE:={}
      start:=i
      repeat
       ROUTE+={start}
       start:=NEXTRETAIL(start)
      until start=i
      ALLRETAIL+=ROUTE
      if getsize(ROUTE)>size then
       LARGEST:=ROUTE
       size:=getsize(LARGEST)
      end-if
     end-if 
    end-do        
   writeln("Selected tour:", LARGEST) 
 
! constraint to exclude the subtour defined by LARGEST
sum(i in LARGEST) ROAD(i,NEXTRETAIL(i)) <= getsize(LARGEST) - 1
 
! force connection from subtour node set to a node not contained in subtour
sum(i in LARGEST,j in RETAIL-LARGEST) ROAD(i,j) >=1
 
! solve the objective function again
   minimize(TOTALDIST)
 
   remove_subtours
  end-if 
 end-procedure