anchor.diff:  !-------------------------------------------------------------------------------
   
   !===============================================================================
  -function st_icosahedron() result(t)
  +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)
  +   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
  
positive.code:    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
     ! 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
  
negative.code:    if (allocated(t%t)) deallocate(t%t)
     allocate(t%t(t%nt))
  
     pole_in = .false.
     if (present(pole)) pole_in = pole
  
     ! First orientation is the division of the octahedron into segments with the golden
     ! ratio, phi
     if (.not.pole) then
        ! 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)
  
        ! 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
  
        ! Normalise onto unit sphere
        call st_norm_p(t%p)
  
     ! Second orientation places two points at +/- z, then places points around small
     ! circles at z = +/-a = 1/sqrt(5).
     ! Coordinates are from http://mathworld.wolfram.com/Icosahedron.html
     else
        t%p(1)  = point(           0._rs,            0._rs,  1._rs)
        t%p(2)  = point(      b*cos(pi5),       b*sin(pi5),      a)
        t%p(3)  = point(b*cos(3._rs*pi5), b*sin(3._rs*pi5),      a)
        t%p(4)  = point(b*cos(5._rs*pi5), b*sin(5._rs*pi5),      a)
        t%p(5)  = point(b*cos(7._rs*pi5), b*sin(7._rs*pi5),      a)
        t%p(6)  = point(b*cos(9._rs*pi5), b*sin(9._rs*pi5),      a)
        t%p(7)  = point(               b,            0._rs,     -a)
        t%p(8)  = point(b*cos(2._rs*pi5), b*sin(2._rs*pi5),     -a)
        t%p(9)  = point(b*cos(4._rs*pi5), b*sin(4._rs*pi5),     -a)
        t%p(10) = point(b*cos(6._rs*pi5), b*sin(6._rs*pi5),     -a)
        t%p(11) = point(b*cos(8._rs*pi5), b*sin(8._rs*pi5),     -a)
        t%p(12) = point(           0._rs,            0._rs, -1._rs)
  
        ! Indices of triangle vertices
        t%t(1)  = triangle( 6,  2,  1)
        t%t(2)  = triangle( 2,  3,  1)
        t%t(3)  = triangle( 3,  4,  1)
        t%t(4)  = triangle( 4,  5,  1)
        t%t(5)  = triangle( 5,  6,  1)
        t%t(6)  = triangle(11,  7,  6)
        t%t(7)  = triangle( 7,  2,  6)
        t%t(8)  = triangle( 7,  8,  2)
        t%t(9)  = triangle( 8,  3,  2)
        t%t(10) = triangle( 8,  9,  3)
        t%t(11) = triangle( 9,  4,  3)
        t%t(12) = triangle( 9, 10,  4)
        t%t(13) = triangle(10,  5,  4)
        t%t(14) = triangle(10, 11,  5)
        t%t(15) = triangle(11,  6,  5)
        t%t(16) = triangle(12,  7, 11)
        t%t(17) = triangle(12,  8,  7)
        t%t(18) = triangle(12,  9,  8)
        t%t(19) = triangle(12, 10,  9)
        t%t(20) = triangle(12, 11, 10)
     endif
  
  end function st_icosahedron
  !-------------------------------------------------------------------------------
  