radiation_force_p Subroutine

  subroutine radiation_force_p(pchunk,q_cache,m_cache,U_X,U_Y,U_Z,E_X,E_Y,E_Z, &
       B_X,B_Y,B_Z,Frad_X,Frad_Y,Frad_Z)
    INTEGER, INTENT(IN)  :: pchunk
    REAL(rp), INTENT(IN)                       :: m_cache,q_cache

    REAL(rp), DIMENSION(pchunk), INTENT(IN)     :: U_X,U_Y,U_Z
    !! \(\mathbf{u} = \gamma \mathbf{v}\), where \(\mathbf{v}\) is the
    !! particle's velocity.
    REAL(rp), DIMENSION(pchunk), INTENT(IN)     :: E_X,E_Y,E_Z
    !! Electric field \(\mathbf{E}\) seen by each particle. This is given
    !! in Cartesian coordinates.
    REAL(rp), DIMENSION(pchunk), INTENT(IN)     :: B_X,B_Y,B_Z
    !! Magnetic field \(\mathbf{B}\) seen by each particle. This is given
    !! in Cartesian coordinates.
    REAL(rp), DIMENSION(pchunk), INTENT(OUT)    :: Frad_X,Frad_Y,Frad_Z
    !! The calculated synchrotron radiation reaction force \(\mathbf{F}_R\).
    REAL(rp), DIMENSION(3)                 :: F1
    !! The component \(\mathbf{F}_1\) of \(\mathbf{F}_R\).
    REAL(rp), DIMENSION(pchunk)                 :: F2_X,F2_Y,F2_Z
    !! The component \(\mathbf{F}_2\) of \(\mathbf{F}_R\).
    REAL(rp), DIMENSION(pchunk)                 :: F3_X,F3_Y,F3_Z
    !! The component \(\mathbf{F}_3\) of \(\mathbf{F}_R\).
    REAL(rp), DIMENSION(pchunk)                 :: V_X,V_Y,V_Z
    !! The particle's velocity \(\mathbf{v}\).
    REAL(rp), DIMENSION(pchunk)                 :: vec_X,vec_Y,vec_Z
    REAL(rp), DIMENSION(pchunk)                 :: cross_EB_X,cross_EB_Y,cross_EB_Z
    REAL(rp), DIMENSION(pchunk)                 :: cross_BV_X,cross_BV_Y,cross_BV_Z
    REAL(rp), DIMENSION(pchunk)                 :: cross_BBV_X,cross_BBV_Y,cross_BBV_Z
    REAL(rp), DIMENSION(pchunk)                 :: dot_EV,dot_vecvec
    !! An auxiliary 3-D vector.
    REAL(rp),DIMENSION(pchunk)                               :: g
    !! The relativistic \(\gamma\) factor of the particle.
    REAL(rp)                               :: tmp
    INTEGER :: cc

    !$OMP SIMD
    !    !$OMP& aligned(g,U_X,U_Y,U_Z,V_X,V_Y,V_Z, &
    !    !$OMP& cross_EB_X,cross_EB_Y,cross_EB_Z,E_X,E_Y,E_Z,B_X,B_Y,B_Z, &
    !    !$OMP& dot_EV,cross_BV_X,cross_BV_Y,cross_BV_Z, &
    !    !$OMP& cross_BBV_X,cross_BBV_Y,cross_BBV_Z,F2_X,F2_Y,F2_Z, &
    !    !$OMP& vec_X,vec_Y,vec_Z,dot_vecvec,F3_X,F3_Y,F3_Z, &
    !    !$OMP& Frad_X,Frad_Y,Frad_Z)
    do cc=1_idef,pchunk
       g(cc) = SQRT(1.0_rp + U_X(cc)*U_X(cc)+ U_Y(cc)*U_Y(cc)+ U_Z(cc)*U_Z(cc))

       V_X(cc) = U_X(cc)/g(cc)
       V_Y(cc) = U_Y(cc)/g(cc)
       V_Z(cc) = U_Z(cc)/g(cc)

       tmp = q_cache**4/(6.0_rp*C_PI*E0*m_cache**2)

       cross_EB_X(cc)=E_Y(cc)*B_Z(cc)-E_Z(cc)*B_Y(cc)
       cross_EB_Y(cc)=E_Z(cc)*B_X(cc)-E_X(cc)*B_Z(cc)
       cross_EB_Z(cc)=E_X(cc)*B_Y(cc)-E_Y(cc)*B_X(cc)

       dot_EV(cc)=E_X(cc)*V_X(cc)+E_Y(cc)*V_Y(cc)+E_Z(cc)*V_Z(cc)

       cross_BV_X(cc)=B_Y(cc)*V_Z(cc)-B_Z(cc)*V_Y(cc)
       cross_BV_Y(cc)=B_Z(cc)*V_X(cc)-B_X(cc)*V_Z(cc)
       cross_BV_Z(cc)=B_X(cc)*V_Y(cc)-B_Y(cc)*V_X(cc)

       cross_BBV_X(cc)=B_Y(cc)*cross_BV_Z(cc)-B_Z(cc)*cross_BV_Y(cc)
       cross_BBV_Y(cc)=B_Z(cc)*cross_BV_X(cc)-B_X(cc)*cross_BV_Z(cc)
       cross_BBV_Z(cc)=B_X(cc)*cross_BV_Y(cc)-B_Y(cc)*cross_BV_X(cc)

       F2_X(cc) = tmp*( dot_EV(cc)*E_X(cc) + cross_EB_X(cc) + cross_BBV_X(cc) )
       F2_Y(cc) = tmp*( dot_EV(cc)*E_Y(cc) + cross_EB_Y(cc) + cross_BBV_Y(cc) )
       F2_Z(cc) = tmp*( dot_EV(cc)*E_Z(cc) + cross_EB_Z(cc) + cross_BBV_Z(cc) )

       vec_X(cc) = E_X(cc) - cross_BV_X(cc)
       vec_Y(cc) = E_Y(cc) - cross_BV_Y(cc)
       vec_Z(cc) = E_Z(cc) - cross_BV_Z(cc)

       dot_vecvec(cc)=vec_X(cc)*vec_X(cc)+vec_Y(cc)*vec_Y(cc)+vec_Z(cc)*vec_Z(cc)

       F3_X(cc) = (tmp*g(cc)**2)*( dot_EV(cc)**2 - dot_vecvec(cc) )*V_X(cc)
       F3_Y(cc) = (tmp*g(cc)**2)*( dot_EV(cc)**2 - dot_vecvec(cc) )*V_Y(cc)
       F3_Z(cc) = (tmp*g(cc)**2)*( dot_EV(cc)**2 - dot_vecvec(cc) )*V_Z(cc)

       Frad_X(cc) = F2_X(cc) + F3_X(cc)
       Frad_Y(cc) = F2_Y(cc) + F3_Y(cc)
       Frad_Z(cc) = F2_Z(cc) + F3_Z(cc)

    end do
    !$OMP END SIMD

  end subroutine radiation_force_p