subroutine calculate_SC_p_FS(params,F,B_R,B_PHI,B_Z,PSIp, &
V_PLL,V_MU,m_cache,flagCon,flagCol,dintJphidPSIP)
TYPE(FIELDS), INTENT(IN) :: F
TYPE(KORC_PARAMS), INTENT(IN) :: params
real(rp),dimension(params%pchunk),intent(in) :: PSIp
real(rp),dimension(params%pchunk),intent(in) :: B_R,B_PHI,B_Z
real(rp),dimension(params%pchunk),intent(in) :: V_PLL,V_MU
real(rp),intent(in) :: m_cache
integer(is),dimension(params%pchunk),intent(in) :: flagCon,flagCol
real(rp),dimension(params%pchunk) :: Bmag,gam,vpll,PSIp_cache
real(rp),dimension(F%dim_1D),intent(out) :: dintJphidPSIP
real(rp),dimension(F%dim_1D) :: PSIP_1D
real(rp) :: dPSIP,ar,arg,arg1,arg2,arg3,PSIP_lim,sigPSIP
integer :: cc,ii,PSIPind,pchunk
pchunk=params%pchunk
PSIP_1D=F%PSIP_1D
dPSIP=PSIP_1D(2)-PSIP_1D(1)
PSIp_cache=PSIp*(params%cpp%Bo*params%cpp%length**2)
sigPSIP=dPSIP
dintJphidPSIP=0._rp
do cc=1_idef,pchunk
! 1D Riemann sum
!write (output_unit_write,*) 'rm',rm(cc)
Bmag(cc)=sqrt(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+ &
B_Z(cc)*B_Z(cc))
gam(cc)=sqrt(1+V_PLL(cc)**2+ &
2*V_MU(cc)*Bmag(cc)*m_cache)
vpll(cc)=V_PLL(cc)/gam(cc)
! write(output_unit_write,*) PSIp_cache(cc)
if (PSIp_cache(cc).lt.0._rp) PSIp_cache(cc)=0._rp
PSIPind=FLOOR(PSIp_cache(cc)/dPSIP)+1_ip
! NGP weighting
! dintJphidPSIP(PSIPind)=dintJphidPSIP(PSIPind)+vpll(cc)
! First-order weighting
! dintJphidPSIP(PSIPind)=dintJphidPSIP(PSIPind)+ &
! vpll(cc)*(PSIP_1D(PSIPind+1)-PSIP_cache(cc))/dPSIP
! dintJphidPSIP(PSIPind+1)=dintJphidPSIP(PSIPind+1)+ &
! vpll(cc)*(PSIP_cache(cc)-PSIP_1D(PSIPind))/dPSIP
! Gaussian weighting
do ii=1_idef,F%dim_1D
arg=MIN((PSIP_1D(ii)-PSIP_cache(cc))**2._rp/ &
(2._rp*sigPSIP**2._rp),100._rp)
dintJphidPSIP(ii)=dintJphidPSIP(ii)+ &
exp(-arg)*vpll(cc)*real(flagCon(cc))*real(flagCol(cc))
end do
end do
! First-order weighting
! dintJphidPSIP(1)=2*dintJphidPSIP(1)
! Gaussian weighting
PSIP_lim=PSIP_1D(F%dim_1D)
do ii=1_idef,F%dim_1D
arg=MIN((PSIP_lim-PSIP_1D(ii))/(sqrt(2._rp)*sigPSIP),10._rp)
arg1=MIN(PSIP_1D(ii)/(sqrt(2._rp)*sigPSIP),10._rp)
dintJphidPSIP(ii)=dintJphidPSIP(ii)/ &
(erf(arg)-erf(-arg1))
end do
end subroutine calculate_SC_p_FS