rpofrho Function

  real function rpofrho(rho)
    implicit none
    real, intent(in) :: rho
    real, parameter :: xerrbi = 1.e-4, xerrsec = 1.e-8
    interface
       real function f(x_)
         implicit none         
         real, intent(in) :: x_
       end function f
    end interface
    procedure(f), pointer :: func
    real :: a, b, fval
    integer :: ier

    ier = 123456 ! Default to indicating error
    a = rpmin ; b = rpmax

    if (irho == 1) then
       func => phi ; fval = rho * rho * func(b)
    else if (irho == 2) then
       func => diameter ; fval = rho * func(b)
    else if (irho == 3) then
       func => psifun ; fval = rho * func(b)
    else if (irho == 4) then
       func => rhofun ; fval = rho * func(b)
    end if

    call root(func, fval, a, b, xerrbi, xerrsec, ier, rpofrho)

    if (debug) then
       write (*,*) "Values in root: "
       write (*,*) 'rho:  ', rho
       write (*,*) 'fval: ', fval
       write (*,*) 'rpofrho: ', rpofrho
       write (*,*) 'rpmin: ', rpmin
       write (*,*) 'rpmax: ', rpmax
    end if
    if (ier > 0) write(*, *) 'error in rpofrho, rho=', rho, rpofrho, ier
  end function rpofrho