anchor.diff:     allocate(t(nt))
      allocate(p(np))
      ! Fill up the new tesselation with the existing np points and nt faces
  -   t(1:u%nt) = u%t(1:u%nt)
  -   p(1:u%np) = u%p(1:u%np)
  -   do it=1,u%nt
  -      t(it)%a = u%t(it)%a;  t(it)%b = u%t(it)%b;  t(it)%c = u%t(it)%c
  -   enddo
  -   do it=1,u%np
  -      p(it)%x = u%p(it)%x;  p(it)%y = u%p(it)%y;  p(it)%z = u%p(it)%z
  -   enddo
  +   t(1:u%nt) = u%t
  +   p(1:u%np) = u%p
      ! Go through each triangle; in each case, we create a new point in the middle
      ! of each of the edges and add these points to the end of the list.
      ! Then replace the original triangle with the one closest to the original point
  
positive.code:    if (allocated(t%t)) deallocate(t%t)
     allocate(t%t(t%nt))
  
     ! Vertices of triangles, unnormalised
     t%p(1)  = point( 0._rs,    phi,  1._rs)
     t%p(2)  = point( 0._rs,   -phi,  1._rs)
     t%p(3)  = point( 0._rs,    phi, -1._rs)
     t%p(4)  = point( 0._rs,   -phi, -1._rs)
     t%p(5)  = point( 1._rs,  0._rs,    phi)
     t%p(6)  = point(-1._rs,  0._rs,    phi)
     t%p(7)  = point( 1._rs,  0._rs,   -phi)
     t%p(8)  = point(-1._rs,  0._rs,   -phi)
     t%p(9)  = point(   phi,  1._rs,  0._rs)
     t%p(10) = point(  -phi,  1._rs,  0._rs)
     t%p(11) = point(   phi, -1._rs,  0._rs)
     t%p(12) = point(  -phi, -1._rs,  0._rs)
  
     ! Normalise onto unit sphere
     call st_norm_p(t%p)
  
     ! Indices describing location of triangle vertices
     t%t( 1)%a =  2 ;  t%t( 1)%b =  4 ;  t%t( 1)%c = 11
     t%t( 2)%a =  5 ;  t%t( 2)%b =  2 ;  t%t( 2)%c = 11
     t%t( 3)%a =  9 ;  t%t( 3)%b =  5 ;  t%t( 3)%c = 11
     t%t( 4)%a =  7 ;  t%t( 4)%b =  9 ;  t%t( 4)%c = 11
     t%t( 5)%a = 11 ;  t%t( 5)%b =  7 ;  t%t( 5)%c =  4
     t%t( 6)%a =  4 ;  t%t( 6)%b =  2 ;  t%t( 6)%c = 12
     t%t( 7)%a =  6 ;  t%t( 7)%b = 12 ;  t%t( 7)%c =  2
     t%t( 8)%a =  2 ;  t%t( 8)%b =  5 ;  t%t( 8)%c =  6
     t%t( 9)%a =  1 ;  t%t( 9)%b =  6 ;  t%t( 9)%c =  5
     t%t(10)%a =  5 ;  t%t(10)%b =  9 ;  t%t(10)%c =  1
     t%t(11)%a =  3 ;  t%t(11)%b =  1 ;  t%t(11)%c =  9
     t%t(12)%a =  9 ;  t%t(12)%b =  7 ;  t%t(12)%c =  3
     t%t(13)%a =  8 ;  t%t(13)%b =  3 ;  t%t(13)%c =  7
     t%t(14)%a =  7 ;  t%t(14)%b =  4 ;  t%t(14)%c =  8
     t%t(15)%a = 12 ;  t%t(15)%b =  8 ;  t%t(15)%c =  4
     t%t(16)%a = 12 ;  t%t(16)%b =  6 ;  t%t(16)%c = 10
     t%t(17)%a =  6 ;  t%t(17)%b =  1 ;  t%t(17)%c = 10
     t%t(18)%a =  1 ;  t%t(18)%b =  3 ;  t%t(18)%c = 10
     t%t(19)%a =  3 ;  t%t(19)%b =  8 ;  t%t(19)%c = 10
     t%t(20)%a = 10 ;  t%t(20)%b = 12 ;  t%t(20)%c =  8
  
  end function st_icosahedron
  !-------------------------------------------------------------------------------
  
negative.code: !-------------------------------------------------------------------------------
  
  !===============================================================================
  function st_icosahedron(pole) result(t)
  !===============================================================================
  !  Return the coordinates and triangles of an icosahedron (12-sided).
  !  By default, the original icosahedron is arranged so that no point lies on the
  !  pole.  This is optimal when looking at data at the poles (see Teanby, C&G, 2006).
  !  Pass pole=.true. to instead use an icosahedron where the poles have points;
  !  this is the arrangement more common in global simulations (see Sadourny,
  !  Arakawa and Mintz, Monthly Weather Review, 1968).
     logical, intent(in), optional :: pole
     logical :: pole_in
     type(tesselation) :: t
     real(rs), parameter :: phi = 2._rs*cos(pi/5._rs), &
        a = 1._rs/sqrt(5._rs), & ! Height of points for polar orientation
        b = 2._rs/sqrt(5._rs), & ! Radius of points for polar orientation
        pi5 = pi/5._rs
  
     ! Number of points
     t%level = 0
  