anchor.diff:      type (meshSphereEquator), intent (in) :: meshEquator
       double precision :: calcDphi
   
  -    if (meshEquator%nPhi .eq. 1) return 0
  -    calcDphi = (meshEquator%phiMax - meshEquator%phiMin) &
  +    if (meshEquator%nPhi .eq. 1) then
  +      calcDphi =  0
  +    else
  +      calcDphi = (meshEquator%phiMax - meshEquator%phiMin) &
                / (meshEquator%nPhi - 1)
  +    end if
     end function calcDphi
   
  -  function calcR(i, meshEquator)
  -    implicit none
  -    integer, intent (in) :: i
  -    type (meshSphereEquator), intent (in) :: meshEquator
  -    double precision :: calcR
  +  function calcR_2(i, meshEquator)
  +    implicit none
  +    integer, intent (in) :: i
  +    type (meshSphereEquator), intent (in) :: meshEquator
  +    double precision :: calcR_2
   
       double precision :: dr
   
       dr = calcDr(meshEquator)
  -    calcR = (i - 0.5) * dr + meshEquator%rMin
  -  end function calcR
  +    calcR_2 = (i - 0.5) * dr + meshEquator%rMin
  +  end function calcR_2
   
     ! A fast version. 
  -  function calcR(i, dr, meshEquator)
  +  function calcR_3(i, dr, meshEquator)
       implicit none
       integer, intent (in) :: i
       double precision, intent (in) :: dr
       type (meshSphereEquator), intent (in) :: meshEquator
  -    double precision :: calcR
  -
  -    calcR = (i - 0.5) * dr + meshEquator%rMin
  -  end function calcR
  -
  -  function calcPhi(i, meshEquator)
  -    implicit none
  -    integer, intent (in) :: i
  -    type (meshSphereEquator), intent (in) :: meshEquator
  -    double precision :: calcPhi
  +    double precision :: calcR_3
  +
  +    calcR_3 = (i - 0.5) * dr + meshEquator%rMin
  +  end function calcR_3
  +
  +  function calcPhi_2(i, meshEquator)
  +    implicit none
  +    integer, intent (in) :: i
  +    type (meshSphereEquator), intent (in) :: meshEquator
  +    double precision :: calcPhi_2
   
       double precision :: dPhi
   
       dPhi = calcDphi(meshEquator)
  -    calcPhi = (i - 0.5) * dPhi + meshEquator%phiMin
  -  end function calcPhi
  +    calcPhi_2 = (i - 0.5) * dPhi + meshEquator%phiMin
  +  end function calcPhi_2
   
     ! A fast version. 
  -  function calcPhi(i, dPhi, meshEquator)
  +  function calcPhi_3(i, dPhi, meshEquator)
       implicit none
       integer, intent (in) :: i
       double precision, intent (in) :: dPhi
       type (meshSphereEquator), intent (in) :: meshEquator
  -    double precision :: calcPhi
  -
  -    calcPhi = (i - 0.5) * dPhi + meshEquator%phiMin
  -  end function calcPhi
  +    double precision :: calcPhi_3
  +
  +    calcPhi_3 = (i - 0.5) * dPhi + meshEquator%phiMin
  +  end function calcPhi_3
   
     ! This is for 2.5D case. 
     subroutine doGrad(res, vals, meshEquator)
  
positive.code:     double precision, intent (in), dimension (:, :) :: Ar, Aphi
      type (meshSphereEquator), intent (in) :: meshEquator
      double precision, intent (out), &
        dimension (size(vals, 1) - 1, size(vals, 2) - 1) :: res
  
      integer :: iR, iPhi, nR, nPhi, nR2, nPhi2
      double precision :: r, dr, dPhi
  
negative.code:     nR = size(vals, 1)
      nPhi = size(vals, 2)
  
      call checkGrid(meshEquator%nR, nR, 'r mesh')
      call checkGrid(meshEquator%nPhi, nPhi, 'phi mesh')
  
      dr = calcDr(meshEquator)
      dPhi = calcDphi(meshEquator)
  