subroutine gravpot(r,nmax,re,mu,c,s,fg)

    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     !! C coefficients
    real(wp),dimension(nmax,0:nmax),intent(in) :: s     !! S coefficients
    real(wp),dimension(3),intent(out)          :: fg    !! 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,temp

    !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 (see p.18)
    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

    main_loop: do i=2,nax0

        if (i/=nax0) then !JW : to prevent access of c,s outside bounds
            ci1 = c(i,1)
            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 in original paper
            !  (should be GOTO 1, rather than 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 main_loop

    !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
            temp = m*a(n,m)
            fac1 = fac1 + temp      *e(n,m)
            fac2 = fac2 + temp      *f(n,m)
            fac3 = fac3 + a(n,m+1)  *d(n,m)
            fac4 = fac4 + a(n+1,m+1)*d(n,m)
        end do

        temp = rinv*rho(n)
        a1 = a1 + temp*fac1
        a2 = a2 + temp*fac2
        a3 = a3 + temp*fac3
        a4 = a4 + temp*fac4

    end do

    fg(1) = a1 - ess*a4
    fg(2) = a2 - t*a4
    fg(3) = a3 - u*a4

    end subroutine gravpot