anchor.diff:  !   Notes:
   !
   !   1.  The employed normalization is the "orthonomalization convention."
  -!   2.  The integral of PlBar**2 over all space on the sphere is 1.
  -!   3.  The integral of PlBar**2 over (-1,1) is 1/2pi.
  +!   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
  
positive.code:         end if
  
       else if (lmax < 0) then
          print*, "Error --- PlBar_d1"
          print*, "LMAX must be greater than or equal to 0."
          print*, "Input value is ", lmax
          if (present(exitstatus)) then
  
negative.code:         end if
  
      else if (abs(z) > 1.0_dp) then
          print*, "Error --- PlON_d1"
          print*, "ABS(Z) must be less than or equal to 1."
          print*, "Input value is ", z
          if (present(exitstatus)) then
  