subroutine grav(mu,rgr,rbar,nmodel,mmodel,cnm,snm,agr)

    implicit none

    real(wp),dimension(3),intent(in)                :: rgr      !! position vector [body-fixed coordinates]
    real(wp),intent(in)                             :: mu       !! gravitational constant
    real(wp),intent(in)                             :: rbar     !! gravitational scaling radius (generally the equatorial radius)
    integer,intent(in)                              :: nmodel   !! the degree of the gravity model (>=2)
    integer,intent(in)                              :: mmodel   !! the order of the gravity model (>=0, <=nmodel)
    real(wp),dimension(nmodel,0:nmodel),intent(in)  :: cnm      !! C gravity coefficients
    real(wp),dimension(nmodel,0:nmodel),intent(in)  :: snm      !! S gravity coefficients
    real(wp),dimension(3),intent(out)               :: agr      !! gravitational acceleration vector [body-fixed coordinates]

    !local variables:
    real(wp),dimension(nmodel,nmodel) :: pnm,ppnm
    real(wp),dimension(nmodel) :: cm,sm,pn,rb,ppn
    real(wp),dimension(3) :: asph
     real(wp) :: e1,e2,e3,e4,e5,r1,r2,t1,t3,absr,sphi,cphi,tcm,tsm,tsnm,tcnm,tpnm
    integer :: n,nm1,nm2,m

    do n=2,nmodel
        pnm(n-1,n) = zero
    end do

    e1 = rgr(1)**2 + rgr(2)**2
    r2 = e1 + rgr(3)**2
    absr = sqrt(r2)
    r1 = sqrt(e1)
    sphi = rgr(3)/absr
    cphi = r1/absr
    if (r1==zero) then
        sm(1) = zero
        cm(1) = one
    else
        sm(1) = rgr(2)/r1
        cm(1) = rgr(1)/r1
    end if
    rb(1) = rbar/absr
    rb(2) = rb(1)**2
    sm(2) = two*cm(1)*sm(1)
    cm(2) = two*cm(1)**2 - one
    pn(1) = sphi
    pn(2) = (three*sphi**2 - one)/two
    ppn(1) = one
    ppn(2) = three*sphi
    pnm(1,1) = one
    pnm(2,2) = three*cphi
    pnm(2,1) = ppn(2)
    ppnm(1,1) = -sphi
    ppnm(2,2) = -six*sphi*cphi
    ppnm(2,1) = three-six*sphi**2

    if (nmodel>=3) then

        do n = 3, nmodel
            nm1 = n-1
            nm2 = n-2
            rb(n) = rb(nm1)*rb(1)
            sm(n) = two*cm(1)*sm(nm1)-sm(nm2)
            cm(n) = two*cm(1)*cm(nm1)-cm(nm2)
            e1 = 2*n-1
            pn(n) = (e1*sphi*pn(nm1)-nm1*pn(nm2))/n
            ppn(n) = sphi*ppn(nm1)+n*pn(nm1)
            pnm(n,n) = e1*cphi*pnm(nm1,nm1)
            ppnm(n,n) = -n*sphi*pnm(n,n)
        end do

        do n = 3, nmodel
            nm1 = n-1
            e1 = (2*n-1)*sphi
            e2 = -n*sphi
            do m = 1, nm1
                e3 = pnm(nm1,m)
                e4 = n+m
                e5 = (e1*e3-(e4-one)*pnm(n-2,m))/(n-m)
                pnm(n,m) = e5
                ppnm(n,m) = e2*e5 + e4*e3
            end do
        end do

    end if

    asph(1) = -one        ![NOTE: set to zero to only output the harmonic terms]
    asph(3) = zero

    do n = 2,nmodel
        e1 = cnm(n,0)*rb(n)
        asph(1) = asph(1)-(n+1)*e1*pn(n)
        asph(3) = asph(3) + e1*ppn(n)
    end do
    asph(3) = cphi*asph(3)
    t1 = zero
    t3 = zero
    asph(2) = zero

    do n = 2, nmodel
        e1 = zero
        e2 = zero
        e3 = zero
        !do m = 1, n                !original
        do m = 1, min(n,mmodel)     !JW - allow for specifying order !!!!!!
            tsnm = snm(n,m)
            tcnm = cnm(n,m)
            tsm = sm(m)
            tcm = cm(m)
            tpnm = pnm(n,m)
            e4 = tsnm*tsm+tcnm*tcm
            e1 = e1+e4*tpnm
            e2 = e2+m*(tsnm*tcm-tcnm*tsm)*tpnm
            e3 = e3+e4*ppnm(n,m)
        end do
        t1 = t1 + (n+1)*rb(n)*e1
        asph(2) = asph(2) + rb(n)*e2
        t3 = t3+rb(n)*e3
    end do

    e4 = mu/r2
    asph(1) = e4*(asph(1) - cphi*t1)
    asph(2) = e4*asph(2)
    asph(3) = e4*(asph(3) + t3)

    e5 = asph(1)*cphi - asph(3)*sphi

    agr(1) = e5*cm(1) - asph(2)*sm(1)
    agr(2) = e5*sm(1) + asph(2)*cm(1)
    agr(3) = asph(1)*sphi + asph(3)*cphi

    end subroutine grav