subroutine geopot(x,y,z,nmax,mmax,re,ksq,c,s,fx,fy,fz)

    implicit none

    !subroutine arguments:
    real(wp),intent(in)              :: x      !! position vector x-component
    real(wp),intent(in)              :: y      !! position vector y-component
    real(wp),intent(in)              :: z      !! position vector z-component
    integer,intent(in)               :: nmax   !! degree of model
    integer,intent(in)               :: mmax   !! order+1 of model
    real(wp),intent(in)              :: re     !! body radius
    real(wp),intent(in)              :: ksq    !! body GM
    real(wp),dimension(:),intent(in) :: c      !! C coefficients
    real(wp),dimension(:),intent(in) :: s      !! S coefficients
    real(wp),intent(out)             :: fx     !! gravitational acceleration x-component
    real(wp),intent(out)             :: fy     !! gravitational acceleration y-component
    real(wp),intent(out)             :: fz     !! gravitational acceleration z-component

    !local variables:
    real(wp) :: r,ri,reor,reorn,ksqor2,xor,yor,zor,rdedx,rdedy,rdedz,&
                sum1,sum2,sum3,sum4,temp1,temp2,temp3,temp4,fact,dcstld,temp
    integer :: i,j,k,im1,l,jm1,jp1,kk
    real(wp),dimension(0:mmax) :: p0,p1,p2,ctil,stil

    !write(*,'(A,1x,*(e20.5,1x/))') 'c=',c
    !write(*,'(A,1x,*(e20.5,1x/))') 's=',s

    !abbreviations:

    r      = sqrt(x*x + y*y + z*z)
    ri     = one/r
    reor   = re*ri
    reorn  = reor
    ksqor2 = ksq*ri*ri
    zor    = z*ri
    xor    = x*ri
    yor    = y*ri

    !the derivatives of the argument of the legendre polynomial - zor

    rdedz = zor*zor - one
    rdedx = zor*xor
    rdedy = zor*yor

    !initialization:

    k = 0
    do i=1,mmax
        p0(i) = zero
        p1(i) = zero
    end do
    p0(0)   = one
    p1(0)   = zor
    p1(1)   = one
    ctil(0) = one
    stil(0) = zero
    ctil(1) = xor
    stil(1) = yor
    sum1    = zero
    !sum2   = zero       !original
    sum2    = one        !JW : include central body term
    sum3    = zero
    sum4    = zero

    !computation of forces:

    do i = 2,nmax

        reorn = reorn*reor
        fact  = 2*i - 1
        im1   = i-1
        l     = 1

        !recursion formulas for legendre polynomial - p2(0)

        p2(0) = (fact*zor*p1(0)-im1*p0(0))/i
        k     = k + 1
        p2(1) = p0(1)+fact*p1(0)
        temp1 = p2(1)*c(k)
        temp2 = p2(0)*c(k)*(i+1)

        if (i < mmax) then

            !recursive formulas for:
            !    'ctilda' - ctil
            !    'stilda' - stil

            ctil(i) = ctil(1)*ctil(im1) - stil(1)*stil(im1)
            stil(i) = stil(1)*ctil(im1) + ctil(1)*stil(im1)
            temp3   = zero
            temp4   = zero

            do j=1,i

                jm1 = j-1
                jp1 = j+1

                !recursive formula for derivative of legendre polynomial - p2(j)

                p2(jp1) = p0(jp1) + fact*p1(j)
                kk      = k + j
                dcstld  = j*p2(j)
                temp    = (c(kk)*ctil(j)+s(kk)*stil(j))
                temp1   = temp1+p2(jp1)*temp
                temp2   = temp2+(i+jp1)*p2(j)*temp
                temp3   = temp3+dcstld*(c(kk)*ctil(jm1)+s(kk)*stil(jm1))
                temp4   = temp4-dcstld*(c(kk)*stil(jm1)-s(kk)*ctil(jm1))

            end do

            l = i
            sum3 = sum3+reorn*temp3
            sum4 = sum4+reorn*temp4

        end if

        sum1 = sum1+reorn*temp1
        sum2 = sum2+reorn*temp2
        k = k + i

        !shift indices:

        do j = 0, l
            p0(j) = p1(j)
            p1(j) = p2(j)
        end do

    end do

    fx = -ksqor2*(sum1*rdedx + sum2*xor - sum3 )
    fy = -ksqor2*(sum1*rdedy + sum2*yor - sum4 )
    fz = -ksqor2*(sum1*rdedz + sum2*zor        )

    end subroutine geopot