pinesnorm Function

    pure function pinesnorm(mu,req,r_f,cnm,snm,nmax,mmax) result(accel)

    implicit none

    real(wp),intent(in)                 :: mu     !! gravitational constant [km^3/s^2]
    real(wp),intent(in)                 :: req    !! body equatorial radius [km]
    real(wp),dimension(3),intent(in)    :: r_f    !! body-fixed Cartesian position vector [km]
    real(wp),dimension(:,0:),intent(in) :: cnm    !! c coefficients (Normalized)
    real(wp),dimension(:,0:),intent(in) :: snm    !! s coefficients (Normalized)
    integer,intent(in)                  :: nmax   !! desired degree
    integer,intent(in)                  :: mmax   !! desired order
    real(wp),dimension(3)               :: accel  !! body-fixed Cartesian acceleration vector [km/s^2]

    real(wp),dimension(nmax+3,nmax+3) :: anm
    real(wp),dimension(mmax+2) :: rm,im
    real(wp),dimension(3) :: harmonics_acc
    real(wp) :: alpha,alpha_den,alpha_num,beta,beta_den,&
                beta_num,dnm,enm,fnm,g1,g1temp,g2,&
                g2temp,g3,g3temp,g4,g4temp,rho,rhop,rmag,&
                s,sm,t,u
    integer :: m,n,m_a,n_a,m_ri,nmodel

    rmag     = norm2(r_f)
    s        = r_f(1)/rmag
    t        = r_f(2)/rmag
    u        = r_f(3)/rmag
    anm      = 0.0_wp
    anm(1,1) = sqrt(2.0_wp)
    do m = 0,nmax+2
        m_a = m + 1
        if (m /= 0) then ! diagonal recursion
           anm(m_a,m_a) = sqrt(1+(1.0_wp/(2*m)))*anm(m_a-1,m_a-1)
        end if
        if (m /= nmax+2) then ! first off-diagonal recursion
            anm(m_a+1,m_a) = sqrt(real(2*m+3,wp))*u*anm(m_a,m_a)
        end if
        if (m < nmax+1) then ! column recursion
            do n = m+2,nmax+2
                n_a          = n + 1
                alpha_num    = (2*n+1)*(2*n-1)
                alpha_den    = (n-m)*(n+m)
                alpha        = sqrt(alpha_num/alpha_den)
                beta_num     = (2*n+1)*(n-m-1)*(n+m-1)
                beta_den     = (2*n-3)*(n+m)*(n-m)
                beta         = sqrt(beta_num/beta_den)
                anm(n_a,m_a) = alpha*u*anm(n_a-1,m_a) - beta*anm(n_a-2,m_a)
            end do
        end if
    end do
    do n = 0,nmax+2
        n_a = n + 1
        anm(n_a,1) = anm(n_a,1)*sqrt(0.5_wp)
    end do
    rm = 0.0_wp
    im = 0.0_wp
    rm(1) = 0.0_wp
    im(1) = 0.0_wp
    rm(2) = 1.0_wp
    im(2) = 0.0_wp
    do m = 1,mmax
        m_ri = m+2
        rm(m_ri) = s*rm(m_ri-1) - t*im(m_ri-1)
        im(m_ri) = s*im(m_ri-1) + t*rm(m_ri-1)
    end do
    rho  = (mu)/(req*rmag)
    rhop = (req)/(rmag)
    g1 = 0.0_wp
    g2 = 0.0_wp
    g3 = 0.0_wp
    g4 = 0.0_wp
    do n = 0,nmax
        n_a = n+1
        g1temp = 0.0_wp
        g2temp = 0.0_wp
        g3temp = 0.0_wp
        g4temp = 0.0_wp
        sm     = 0.5_wp
        if (n>mmax) then
            nmodel=mmax
        else
            nmodel=n
        end if
        do m = 0,nmodel
            if (n>=2) then
                m_a    = m + 1
                m_ri   = m + 2
                dnm    = cnm(n,m)*rm(m_ri)   + snm(n,m)*im(m_ri)
                enm    = cnm(n,m)*rm(m_ri-1) + snm(n,m)*im(m_ri-1)
                fnm    = snm(n,m)*rm(m_ri-1) - cnm(n,m)*im(m_ri-1)
                alpha  = sqrt(sm*(n-m)*(n+m+1))
                g1temp = g1temp + anm(n_a,m_a)*(m)*enm
                g2temp = g2temp + anm(n_a,m_a)*(m)*fnm
                g3temp = g3temp + alpha*anm(n_a,m_a+1)*dnm
                g4temp = g4temp + ((n+m+1)*anm(n_a,m_a)+alpha*u*anm(n_a,m_a+1))*dnm
                if (m == 0) sm = 1.0_wp
            end if
        end do
        rho = rhop*rho
        g1  = g1 + rho*g1temp
        g2  = g2 + rho*g2temp
        g3  = g3 + rho*g3temp
        g4  = g4 + rho*g4temp
    end do

    ! include central body term here:
    accel = (-mu/rmag**3 * r_f) + [g1 - g4*s, g2 - g4*t, g3 - g4*u]

    end function pinesnorm