interp_FOfields_aorsa_p Subroutine

subroutine interp_FOfields_aorsa_p(time,params,pchunk,F,Y_R,Y_PHI,Y_Z, &
     B_X,B_Y,B_Z,E_X,E_Y,E_Z,PSIp,flag_cache)
  TYPE(KORC_PARAMS), INTENT(IN)                              :: params
  REAL(rp), INTENT(IN)  :: time
  INTEGER, INTENT(IN)  :: pchunk
  TYPE(FIELDS), INTENT(IN)                               :: F
  REAL(rp),DIMENSION(pchunk),INTENT(IN)   :: Y_R,Y_PHI,Y_Z
  REAL(rp),DIMENSION(pchunk),INTENT(OUT)   :: B_X,B_Y,B_Z
  REAL(rp),DIMENSION(pchunk),INTENT(OUT)   :: E_X,E_Y,E_Z
  REAL(rp),DIMENSION(pchunk)   :: B0_R,B0_PHI,B0_Z,B0_X,B0_Y
  REAL(rp),DIMENSION(pchunk)   :: B1_X,B1_Y,B1_Z
  REAL(rp),DIMENSION(pchunk)   :: B1Re_X,B1Re_Y,B1Re_Z
  REAL(rp),DIMENSION(pchunk)   :: B1Im_X,B1Im_Y,B1Im_Z
  REAL(rp),DIMENSION(pchunk)   :: E1Re_X,E1Re_Y,E1Re_Z
  REAL(rp),DIMENSION(pchunk)   :: E1Im_X,E1Im_Y,E1Im_Z
  REAL(rp),DIMENSION(pchunk),INTENT(OUT)   :: PSIp
  REAL(rp),DIMENSION(pchunk)   :: cP,sP,cnP,snP
  REAL(rp), DIMENSION(pchunk,3)  :: A
  !  INTEGER(ip) :: ezerr
  INTEGER                                      :: cc
  !! Particle chunk iterator.
  INTEGER(is),DIMENSION(pchunk),INTENT(INOUT)   :: flag_cache
  REAL(rp) :: psip_conv
  REAL(rp) :: amp,nmode,omega

  psip_conv=F%psip_conv
  amp=F%AMP
  nmode=F%AORSA_nmode
  omega=2*C_PI*F%AORSA_freq
  
  call check_if_in_fields_domain_p(pchunk,F,Y_R,Y_PHI,Y_Z,flag_cache)

  call EZspline_interp(bfield_2d%A,b1Refield_2dx%X,b1Refield_2dx%Y, &
       b1Refield_2dx%Z,b1Imfield_2dx%X,b1Imfield_2dx%Y,b1Imfield_2dx%Z, &
       e1Refield_2dx%X,e1Refield_2dx%Y,e1Refield_2dx%Z, &
       e1Imfield_2dx%X,e1Imfield_2dx%Y,e1Imfield_2dx%Z, &
       pchunk,Y_R,Y_Z,A,B1Re_X,B1Re_Y,B1Re_Z,B1Im_X,B1Im_Y,B1Im_Z, &
       E1Re_X,E1Re_Y,E1Re_Z,E1Im_X,E1Im_Y,E1Im_Z,ezerr)
  call EZspline_error(ezerr)

  
  !$OMP SIMD
  do cc=1_idef,pchunk
     PSIp(cc)=A(cc,1)

     B0_R(cc) = psip_conv*A(cc,3)/Y_R(cc)
     B0_PHI(cc) = -F%Bo*F%Ro/Y_R(cc)
     B0_Z(cc) = -psip_conv*A(cc,2)/Y_R(cc)

     cP(cc)=cos(Y_PHI(cc))
     sP(cc)=sin(Y_PHI(cc))
     
     B0_X(cc) = B0_R(cc)*cP(cc) - B0_PHI(cc)*sP(cc)
     B0_Y(cc) = B0_R(cc)*sP(cc) + B0_PHI(cc)*cP(cc)
     

     cnP(cc)=cos(omega*time+nmode*Y_PHI(cc))
     snP(cc)=sin(omega*time+nmode*Y_PHI(cc))
     
     B1_X(cc) = amp*(B1Re_X(cc)*cnP(cc)-B1Im_X(cc)*snP(cc))
     B1_Y(cc) = amp*(B1Re_Y(cc)*cnP(cc)-B1Im_Y(cc)*snP(cc))
     B1_Z(cc) = amp*(B1Re_Z(cc)*cnP(cc)-B1Im_Z(cc)*snP(cc))

     B_X(cc) = B0_X(cc)+B1_X(cc)
     B_Y(cc) = B0_Y(cc)+B1_Y(cc)
     B_Z(cc) = B0_Z(cc)+B1_Z(cc)
     
     E_X(cc) = amp*(E1Re_X(cc)*cnP(cc)-E1Im_X(cc)*snP(cc))
     E_Y(cc) = amp*(E1Re_Y(cc)*cnP(cc)-E1Im_Y(cc)*snP(cc))
     E_Z(cc) = amp*(E1Re_Z(cc)*cnP(cc)-E1Im_Z(cc)*snP(cc))     


  end do
  !$OMP END SIMD

#if DBG_CHECK    
  !write(6,*) '(R,PHI,Z,time)',Y_R*params%cpp%length,Y_PHI, &
  !     Y_Z*params%cpp%length,time
  !write(6,*) 'psi',PSIp*params%cpp%Bo*params%cpp%length**2
  !write(6,*) 'dpsidR',A(:,2)*params%cpp%Bo*params%cpp%length
  !write(6,*) 'dpsidZ',A(:,3)*params%cpp%Bo*params%cpp%length
  !write(6,*) 'B0',B0_R*params%cpp%Bo,B0_PHI*params%cpp%Bo,B0_Z*params%cpp%Bo
  !write(6,*) 'AMP',amp
  !write(6,*) 'B1Re',B1Re_X*params%cpp%Bo,B1Re_Y*params%cpp%Bo,B1Re_Z*params%cpp%Bo
  !write(6,*) 'B1Im',B1Im_X*params%cpp%Bo,B1Im_Y*params%cpp%Bo,B1Im_Z*params%cpp%Bo
  !write(6,*) 'B1',B1_X*params%cpp%Bo,B1_Y*params%cpp%Bo,B1_Z*params%cpp%Bo
  !write(6,*) 'B',B_X*params%cpp%Bo,B_Y*params%cpp%Bo,B_Z*params%cpp%Bo
#endif
  
end subroutine interp_FOfields_aorsa_p