subroutine calculate_GCfields_p_FS(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z, &
E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, &
flag_cache,PSIp)
INTEGER, INTENT(IN) :: pchunk
REAL(rp), DIMENSION(pchunk), INTENT(IN) :: Y_R,Y_PHI,Y_Z
TYPE(FIELDS), INTENT(IN) :: F
REAL(rp), DIMENSION(pchunk), INTENT(OUT) :: B_R,B_PHI,B_Z
REAL(rp), DIMENSION(pchunk), INTENT(OUT) :: gradB_R,gradB_PHI,gradB_Z
REAL(rp), DIMENSION(pchunk), INTENT(OUT) :: curlb_R,curlb_PHI,curlb_Z
REAL(rp), DIMENSION(pchunk), INTENT(OUT) :: E_R,E_PHI,E_Z
REAL(rp), DIMENSION(pchunk) :: Bmag
INTEGER :: cc
REAL(rp), DIMENSION(pchunk),INTENT(OUT) :: PSIp
REAL(rp), DIMENSION(pchunk,6) :: A
INTEGER(is),DIMENSION(pchunk),INTENT(INOUT) :: flag_cache
REAL(rp) :: psip_conv
psip_conv=F%psip_conv
call check_if_in_fields_domain_p(pchunk,F,Y_R,Y_PHI,Y_Z,flag_cache)
call EZspline_derivative(bfield_2d%A, pchunk, Y_R, Y_Z, A, ezerr)
call EZspline_error(ezerr)
!write (output_unit_write,*) A(1,1),A(1,2)
!A(:,1) = PSIp
!A(:,2) = dPSIp/dR
!A(:,3) = dPSIp/dZ
!A(:,4) = d^2PSIp/dR^2
!A(:,5) = d^2PSIp/dZ^2
!A(:,6) = d^2PSIp/dRdZ
!$OMP SIMD
! !$OMP& aligned(PSIp,A,B_R,Y_R,B_PHI,B_Z,Bmag,gradB_R,gradB_PHI,gradB_Z, &
! !$OMP& curlb_R,curlb_PHI,curlb_Z,E_R,E_PHI,E_Z)
do cc=1_idef,pchunk
PSIp(cc)=A(cc,1)
A(cc,2)=psip_conv*A(cc,2)/(2*C_PI)
A(cc,3)=psip_conv*A(cc,3)/(2*C_PI)
A(cc,4)=psip_conv*A(cc,4)/(2*C_PI)
A(cc,5)=psip_conv*A(cc,5)/(2*C_PI)
A(cc,6)=psip_conv*A(cc,6)/(2*C_PI)
B_R(cc) = A(cc,3)/Y_R(cc)
! BR = (dA/dZ)/R
B_PHI(cc) = -F%Bo*F%Ro/Y_R(cc)
! BPHI = Fo*Ro/R
B_Z(cc) = -A(cc,2)/Y_R(cc)
! BR = -(dA/dR)/R
Bmag(cc)=sqrt(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+B_Z(cc)*B_Z(cc))
gradB_R(cc)=(B_R(cc)*A(cc,6)-B_Z(cc)*A(cc,4)-Bmag(cc)*Bmag(cc))/ &
(Y_R(cc)*Bmag(cc))
gradB_PHI(cc)=0._rp
gradB_Z(cc)=(B_R(cc)*A(cc,5)-B_Z(cc)*A(cc,6))/ &
(Y_R(cc)*Bmag(cc))
curlb_R(cc)=B_PHI(cc)*gradB_Z(cc)/(Bmag(cc)*Bmag(cc))
curlb_PHI(cc)=(Bmag(cc)/Y_R(cc)*(B_Z(cc)+A(cc,4)+A(cc,5))- &
B_R(cc)*gradB_Z(cc)+B_Z(cc)*gradB_R(cc))/ &
(Bmag(cc)*Bmag(cc))
curlb_Z(cc)=-B_PHI(cc)*gradB_R(cc)/(Bmag(cc)*Bmag(cc))
E_R(cc) = 0._rp
E_PHI(cc) = F%Eo*F%Ro/Y_R(cc)
E_Z(cc) = 0._rp
end do
end subroutine calculate_GCfields_p_FS