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: 
      if (z == 1.0_dp) then
          do l = 0, lmax
              p(l+1) = sqrt( dble(2*l+1) )
              dp1(l+1) = sqrt( dble(2*l+1) ) * dble(l) * dble(l+1) / 2.0_dp
          end do
  
      else if (z == -1.0_dp) then
          do l = 0, lmax
              p(l+1) = sqrt( dble(2*l+1) ) * dble((-1)**l)
              dp1(l+1) = sqrt( dble(2*l+1)) * dble(l) * dble(l+1) &
                          * dble((-1)**(l-1)) / 2.0_dp
              end do
  
      else
          sinsq = (1.0_dp - z**2)
  
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
  