subroutine calculate_3DBdBfields_p(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)
INTEGER, INTENT(IN) :: pchunk
REAL(rp), DIMENSION(pchunk), INTENT(IN) :: Y_R,Y_PHI,Y_Z
real(rp), DIMENSION(pchunk) :: Y_PHI_mod
TYPE(FIELDS), INTENT(IN) :: F
REAL(rp), DIMENSION(pchunk), INTENT(OUT) :: B_R,B_PHI,B_Z
REAL(rp), DIMENSION(pchunk) :: dBRdR,dBPHIdR,dBZdR
REAL(rp), DIMENSION(pchunk) :: dBRdPHI,dBPHIdPHI,dBZdPHI
REAL(rp), DIMENSION(pchunk) :: dBRdZ,dBPHIdZ,dBZdZ
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
INTEGER(is),DIMENSION(pchunk),INTENT(INOUT) :: flag_cache
Y_PHI_mod=modulo(Y_PHI,2._rp*C_PI)
call check_if_in_fields_domain_p(pchunk,F,Y_R,Y_PHI_mod,Y_Z,flag_cache)
call EZspline_interp(bfield_2d%R,bfield_2d%PHI,bfield_2d%Z, &
dbdR_2d%R,dbdR_2d%PHI,dBdR_2d%Z, &
dbdPHI_2d%R,dbdPHI_2d%PHI,dbdPHI_2d%Z, &
dbdZ_2d%R,dbdZ_2d%PHI,dbdZ_2d%Z, &
efield_2d%R,efield_2d%PHI,efield_2d%Z,pchunk,Y_R,Y_Z,B_R,B_PHI,B_Z, &
dBRdR,dBPHIdR,dBZdR,dBRdPHI,dBPHIdPHI,dBZdPHI,dBRdZ,dBPHIdZ,dBZdZ, &
E_R,E_PHI,E_Z,ezerr)
call EZspline_error(ezerr)
!$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
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)*dBRdR(cc)+B_PHI(cc)*dBPHIdR(cc)+ &
B_Z(cc)*dBZdR(cc))/Bmag(cc)
gradB_PHI(cc)=(B_R(cc)*dBRdPHI(cc)+B_PHI(cc)*dBPHIdPHI(cc)+ &
B_Z(cc)*dBZdPHI(cc))/(Y_R(cc)*Bmag(cc))
gradB_Z(cc)=(B_R(cc)*dBRdZ(cc)+B_PHI(cc)*dBPHIdZ(cc)+ &
B_Z(cc)*dBZdZ(cc))/Bmag(cc)
curlb_R(cc)=(Bmag(cc)*dBZdPHI(cc)/Y_R(cc)-B_Z(cc)*gradB_PHI(cc)- &
Bmag(cc)*dBPHIdZ(cc)+B_PHI(cc)*gradB_Z(cc))/(Bmag(cc)*Bmag(cc))
curlb_PHI(cc)=(Bmag(cc)*dBRdZ(cc)-B_R(cc)*gradB_Z(cc)- &
Bmag(cc)*dBZdR(cc)+B_Z(cc)*gradB_R(cc))/(Bmag(cc)*Bmag(cc))
curlb_Z(cc)=(Bmag(cc)*B_PHI(cc)/Y_R(cc)+Bmag(cc)*dBPHIdR(cc)- &
B_PHI(cc)*gradB_R(cc)-Bmag(cc)*dBRdPHI(cc)/Y_R(cc)+ &
B_R(cc)*gradB_PHI(cc))/(Bmag(cc)*Bmag(cc))
end do
!$OMP END SIMD
! write(output_unit_write,'("PSIp: ",E17.10)') PSIp
! write(output_unit_write,'("Y_R: ",E17.10)') Y_R
! write(output_unit_write,'("Y_Z: ",E17.10)') Y_Z
! write(output_unit_write,'("B_R: ",E17.10)') B_R
! write(output_unit_write,'("B_PHIinterp: ",E17.10)') B_PHI
! write(output_unit_write,'("B_Z: ",E17.10)') B_Z
end subroutine calculate_3DBdBfields_p