$$
U = \frac{\mu}{r} [ 1 + \sum\limits_{n=2}^{n_{max}} \sum\limits_{m=0}^{n} \left( \frac{r_e}{r} \right)^n P_{n,m}(\sin \phi) ( C_{n,m} \cos m\lambda + S_{n,m} \sin m\lambda ) ]
$$
A geopotential formulation that eliminates the polar singularity was devised by Samuel Pines in the early 1970s [1]. Over the years, various implementations and refinements of this algorithm and other similar algorithms have been published [2-7], mostly at the Johnson Space Center. The Fortran 77 code in Mueller [2] and Lear [4-5] is easily translated into modern Fortran. The Pines algorithm code from Spencer [3] appears to contain a few bugs, so translation is more of a challenge. A modern Fortran version (with bugs fixed) is presented below:
subroutine gravpot(r,nmax,re,mu,c,s,acc)
use, intrinsic :: iso_fortran_env, only: wp => real64
implicit none
real(wp),dimension(3),intent(in) :: r !position vector
integer,intent(in) :: nmax !degree/order
real(wp),intent(in) :: re !body radius
real(wp),intent(in) :: mu !grav constant
real(wp),dimension(nmax,0:nmax),intent(in) :: c !coefficients
real(wp),dimension(nmax,0:nmax),intent(in) :: s !
real(wp),dimension(3),intent(out) :: acc !grav acceleration
!local variables:
real(wp),dimension(nmax+1) :: creal, cimag, rho
real(wp),dimension(nmax+1,nmax+1) :: a,d,e,f
integer :: nax0,i,j,k,l,n,m
real(wp) :: rinv,ess,t,u,r0,rhozero,a1,a2,a3,a4,fac1,fac2,fac3,fac4
real(wp) :: ci1,si1
real(wp),parameter :: zero = 0.0_wp
real(wp),parameter :: one = 1.0_wp
!JW : not done in original paper,
! but seems to be necessary
! (probably assumed the compiler
! did it automatically)
a = zero
d = zero
e = zero
f = zero
!get the direction cosines ess, t and u:
nax0 = nmax + 1
rinv = one/norm2(r)
ess = r(1) * rinv
t = r(2) * rinv
u = r(3) * rinv
!generate the functions creal, cimag, a, d, e, f and rho:
r0 = re*rinv
rhozero = mu*rinv !JW: typo in original paper
rho(1) = r0*rhozero
creal(1) = ess
cimag(1) = t
d(1,1) = zero
e(1,1) = zero
f(1,1) = zero
a(1,1) = one
do i=2,nax0
if (i/=nax0) then !JW : to prevent access
ci1 = c(i,1) ! to c,s outside bounds
si1 = s(i,1)
else
ci1 = zero
si1 = zero
end if
rho(i) = r0*rho(i-1)
creal(i) = ess*creal(i-1) - t*cimag(i-1)
cimag(i) = ess*cimag(i-1) + t*creal(i-1)
d(i,1) = ess*ci1 + t*si1
e(i,1) = ci1
f(i,1) = si1
a(i,i) = (2*i-1)*a(i-1,i-1)
a(i,i-1) = u*a(i,i)
do k=2,i
if (i/=nax0) then
d(i,k) = c(i,k)*creal(k) + s(i,k)*cimag(k)
e(i,k) = c(i,k)*creal(k-1) + s(i,k)*cimag(k-1)
f(i,k) = s(i,k)*creal(k-1) - c(i,k)*cimag(k-1)
end if
!JW : typo here in original paper
! (should be GOTO 1, not GOTO 10)
if (i/=2) then
L = i-2
do j=1,L
a(i,i-j-1) = (u*a(i,i-j)-a(i-1,i-j))/(j+1)
end do
end if
end do
end do
!compute auxiliary quantities a1, a2, a3, a4
a1 = zero
a2 = zero
a3 = zero
a4 = rhozero*rinv
do n=2,nmax
fac1 = zero
fac2 = zero
fac3 = a(n,1) *c(n,0)
fac4 = a(n+1,1)*c(n,0)
do m=1,n
fac1 = fac1 + m*a(n,m) *e(n,m)
fac2 = fac2 + m*a(n,m) *f(n,m)
fac3 = fac3 + a(n,m+1) *d(n,m)
fac4 = fac4 + a(n+1,m+1)*d(n,m)
end do
a1 = a1 + rinv*rho(n)*fac1
a2 = a2 + rinv*rho(n)*fac2
a3 = a3 + rinv*rho(n)*fac3
a4 = a4 + rinv*rho(n)*fac4
end do
!gravitational acceleration vector:
acc(1) = a1 - ess*a4
acc(2) = a2 - t*a4
acc(3) = a3 - u*a4
end subroutine gravpot
