subroutine kuga_carrara_geopotential(nmax, nm, re, gm, c, s, x, ac)

    implicit none

    integer,intent(in)                       :: nmax !! max. order and degree loaded
    integer,intent(in)                       :: nm   !! desired order and degree (nm <= nmax)
    real(wp),intent(in)                      :: re   !! body equatorial radius [km]
    real(wp),intent(in)                      :: gm   !! gravitational constant [km3/s2]
    real(wp),dimension(nmax,0:nm),intent(in) :: c    !! c coefficients (Normalized)
    real(wp),dimension(nmax,0:nm),intent(in) :: s    !! s coefficients (Normalized)
    real(wp),dimension(3),intent(in)         :: x    !! body-fixed cartesian position vector [km]
    real(wp),dimension(3),intent(out)        :: ac   !! body-fixed cartesian acceleration vector [km/s2]

    integer :: n, m
    real(wp),dimension(0:nm) :: pn, qn
    real(wp) :: r, q, t, u2, u, um, tf, al, sl, cl, gmr
    real(wp) :: pnm, dpnm, anm, bnm, fnm, cmm, smm
    real(wp) :: am, an, pnn, pnm1m, pnm2m, sm, cm, sml, cml
    real(wp) :: qc, qs, xc, xs, xcf, xsf, xcr, xsr, vl, vf, vr

    real(wp),parameter :: sqrt3 = sqrt(3.0_wp)

    ! auxiliary variables
    r   = norm2(x)
    q   = re / r
    t   = x(3) / r             ! sin (lat)
    u   = sqrt(1.0_wp - t*t)
    tf  = t/u                  ! tan (lat)
    al  = atan2(x(2), x(1))
    sl  = sin (al)             ! sin (long)
    cl  = cos (al)             ! cos (long)
    gmr = gm / r

    ! initialize
    vl = 0.0_wp
    vf = 0.0_wp
    vr = 0.0_wp

    ! store sectoral
    pn(0) = 1.0_wp
    pn(1) = sqrt3 * u ! sqrt(3) * cos (lat)
    qn(0) = 1.0_wp
    qn(1) = q
    do m = 2, nm
        am = real(m,wp)
        pn(m) = u * sqrt(1.0_wp+0.5_wp/am) * pn(m-1)
        qn(m) = q * qn(m-1)
    end do

    ! initialize sin and cos recursions
    sm = 0.0_wp
    cm = 1.0_wp

    ! outer loop
    do m = 0,nm

        ! init
        am = real(m,wp)

        ! for m = n (sectoral)
        pnm   = pn(m)
        dpnm  = -am * pnm * tf
        pnm1m = pnm
        pnm2m = 0.0_wp

        ! initialize Horner's scheme
        if (m<2) then
            cmm = 0.0_wp
            smm = 0.0_wp
        else
            cmm = c(m,m)
            smm = s(m,m)
        end if

        qc  = qn(m) * cmm
        qs  = qn(m) * smm
        xc  = qc * pnm
        xs  = qs * pnm
        xcf = qc * dpnm
        xsf = qs * dpnm
        xcr = (am+1.0_wp) * qc * pnm
        xsr = (am+1.0_wp) * qs * pnm

        ! inner loop
        do n = m+1,nm

            an  = real(n,wp)
            anm = sqrt(((an+an-1.0_wp)*(an+an+1.0_wp))/((an-am)*(an+am)))
            bnm = sqrt(((an+an+1.0_wp)*(an+am-1.0_wp)*(an-am-1.0_wp))/((an-am)*(an+am)*(an+an-3.0_wp)))
            fnm = sqrt(((an*an-am*am)*(an+an+1.0_wp))/(an+an-1.0_wp))

            ! recursion p and dp
            pnm  = anm * t * pnm1m - bnm * pnm2m
            dpnm = -an * pnm * tf + fnm * pnm1m / u   ! signal opposite to paper

            ! store
            pnm2m = pnm1m
            pnm1m = pnm

            ! inner sum
            if (n >= 2) then
                qc  = qn(n) * c(n,m)
                qs  = qn(n) * s(n,m)
                xc  = xc + qc * pnm
                xs  = xs + qs * pnm
                xcf = xcf + qc * dpnm
                xsf = xsf + qs * dpnm
                xcr = xcr + (an+1.0_wp) * qc * pnm
                xsr = xsr + (an+1.0_wp) * qs * pnm
            end if

        end do

        ! outer sum
        vl = vl + am*(xc*sm - xs*cm)
        vf = vf + (xcf*cm + xsf*sm)
        vr = vr + (xcr*cm + xsr*sm)

        ! sin and cos recursions for next m
        cml = cl*cm - sm*sl
        sml = cl*sm + cm*sl
        cm  = cml ! save for next m
        sm  = sml ! save for next m

    end do

    ! gradient
    vl = -gmr * vl
    vf =  gmr * vf
    vr = -(gmr/r) * (1.0_wp+vr)

    ! body x, y, z accelerations
    ac(1) = u*cl*vr - t*cl*vf/r - sl*vl/(u*r)
    ac(2) = u*sl*vr - t*sl*vf/r + cl*vl/(u*r)
    ac(3) = t*vr + u*vf/r

    end subroutine kuga_carrara_geopotential