anchor.diff:          end if
   
        else if (lmax < 0) then
  -        print*, "Error --- PlBar_d1"
  +        print*, "Error --- PlON_d1"
           print*, "LMAX must be greater than or equal to 0."
           print*, "Input value is ", lmax
           if (present(exitstatus)) then
  
positive.code:     pi = acos(-1.0_dp)
  
      if (size(p) < lmax+1) then
          print*, "Error --- PlBar_d1"
          print*, "P must be dimensioned as (LMAX+1) where LMAX is ", lmax
          print*, "Input array is dimensioned ", size(p)
          if (present(exitstatus)) then
  
negative.code: !   Notes:
  !
  !   1.  The employed normalization is the "orthonomalization convention."
  !   2.  The integral of PlON**2 over all space on the sphere is 1.
  !   3.  The integral of PlON**2 over (-1,1) is 1/2pi.
  !   4.  The derivative is evaluated with respecte to z, and NOT cos(colatitude)
  !       or sin(latitude).
  !   5.  Derivatives are calculated according to the geodesy-normalized
  