calculate_magnetic_field_p Subroutine

subroutine calculate_magnetic_field_p(F,Y_R,Y_Z,B_R,B_PHI,B_Z)
  REAL(rp), DIMENSION(8), INTENT(IN)      :: Y_R,Y_Z
  TYPE(FIELDS), INTENT(IN)                               :: F
  REAL(rp), DIMENSION(8),  INTENT(OUT)   :: B_R,B_PHI,B_Z
  INTEGER                                                :: pp
  REAL(rp), DIMENSION(8)  :: PSIp
  REAL(rp), DIMENSION(8,2)  :: A
  
  call EZspline_interp(bfield_2d%A, 8, Y_R, Y_Z, &
       PSIp, ezerr)
  call EZspline_error(ezerr)
     
  ! FR = (dA/dZ)/R
  call EZspline_gradient(bfield_2d%A, 8, Y_R, Y_Z, &
       A, ezerr)
  call EZspline_error(ezerr)

  !write(6,'("dPSIp/dR: ",E17.10)') A(:,1)
  !write(6,'("dPSIp/dZ: ",E17.10)') A(:,2)
  !write(6,'("Y_R: ",E17.10)') Y_R

  B_R = A(:,2)/Y_R

  ! FPHI = Fo*Ro/R
  B_PHI = -F%Bo*F%Ro/Y_R

  ! FR = -(dA/dR)/R

  !     write(6,'("R*B_Z: ",E17.10)') B_Z(1)

  B_Z= -A(:,1)/Y_R

   
!  write(6,'("PSIp: ",E17.10)') PSIp
  
!  write(6,'("Y_R: ",E17.10)') Y_R(1)
!  write(6,'("Y_Z: ",E17.10)') Y_Z(1)
  
!  write(6,'("B_R: ",E17.10)') B_R
!  write(6,'("B_PHI: ",E17.10)') B_PHI
!  write(6,'("B_Z: ",E17.10)') B_Z

end subroutine calculate_magnetic_field_p