anchor.diff: +!===============================================================================
  +! One of Andy Nowacki's Fortran utility modules for dealing with seismic
  +! anisotropy and other problems.
  +!
  +! Andy Nowacki <andy.nowacki@bristol.ac.uk>
  +!
  +! See the file LICENCE for licence details.
   !===============================================================================
   program sph_tess
   !===============================================================================
   ! Dump the points for a level-n tesselation of a sphere, starting with an
   ! icosahedron.
  +! This version will create a tesselation of higher order if it is not available
  +! as a pre-saved version on disk.
   
      use sphere_tesselate
      use spherical_geometry
  
positive.code:    implicit none
     type(tesselation) :: t
     integer :: i, n
     logical :: output_triangles = .false., output_geog = .false.
  
     call get_args
  
     t = st_icosahedron(pole=.true.)
     do i = 1, n
        call st_iterate_level(t)
     enddo
     if (.not.output_triangles) then
        call st_dump_points(t, geog=output_geog)
     else
  
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
  
  