anchor.diff:     implicit none
      type(tesselation) :: t
      integer :: i, n
  -   logical :: output_triangles = .false., output_geog = .false.
  +   logical :: output_triangles = .false., output_geog = .false., &
  +      cache_exists = .false.
  +   ! Path to directory containing pre-computed tesselations
  +   character(len=250) :: file
   
      call get_args
   
  -   t = st_icosahedron(pole=.true.)
  -   do i = 1, n
  -      call st_iterate_level(t)
  -   enddo
  +   ! Try to load the file from the global cache
  +   call st_load_cache(t, cache_exists)
  +   ! Otherwise, create the tesselation and save the cache file for later use
  +   if (.not.cache_exists) then
  +      write(0,'(a)') 'Generating tesselation and saving cache file...'
  +      t = st_icosahedron(pole=.true.)
  +      do i = 1, n
  +         call st_iterate_level(t)
  +      enddo
  +      call st_save_cache(t)
  +   endif
  +   ! Output
      if (.not.output_triangles) then
         call st_dump_points(t, geog=output_geog)
      else
  
positive.code:          'icosahedron.  x, y, z values on the unit sphere are output to stdout.', &
           'n must be >= 0', &
           'Options:', &
           '   -g : Write values as lon, lat instead of x, y, z', &
           '   -t : Dump the triangles as multisegment lines'
        stop
     end subroutine usage
  
  
negative.code:       do while (iarg <= narg - 1)
           call get_command_argument(iarg, arg)
           select case(arg)
              case('-g');   output_geog = .true.;       iarg = iarg + 1
              case('-p');   call get_command_argument(iarg+1, st_cache_dir)
                                                        iarg = iarg + 2
              case('-t');   output_triangles = .true.;  iarg = iarg + 1
              case default; call usage
           end select
        enddo
        ! Get tesselation level
        call get_command_argument(narg, arg)
        read(arg,*) n
        t%level = n
        if (n < 0) call usage
     end subroutine get_args
  
  