calculate_SC_p Subroutine

  subroutine calculate_SC_p(params,F,B_R,B_PHI,B_Z,Y_R,Y_Z, &
       V_PLL,V_MU,m_cache,flagCon,flagCol,Vden)

    TYPE(FIELDS), INTENT(IN)                 :: F
    TYPE(KORC_PARAMS), INTENT(IN) 		:: params
    real(rp),dimension(params%pchunk),intent(in) :: Y_R,Y_Z
    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) :: rm,Bmag,gam,vpll
    real(rp),dimension(F%dim_1D),intent(out) :: Vden
    real(rp),dimension(F%dim_1D) :: Vpart,Ai
    real(rp),dimension(F%dim_1D) :: r_1D
    real(rp) :: dr,sigr,ar,arg,arg1,arg2,arg3
    integer :: cc,ii,rind,pchunk

    pchunk=params%pchunk

    dr=F%r_1D(2)-F%r_1D(1)
    r_1D=F%r_1D
    sigr=dr

    Vpart=0._rp
    do cc=1_idef,pchunk

       ! 1D nearest grid point weighting in minor radius

       !    RR=spp%vars%Y(:,1)
       !    ZZ=spp%vars%Y(:,3)
       rm(cc)=sqrt((Y_R(cc)-F%Ro)**2+(Y_Z(cc)-F%Zo)**2)* &
            params%cpp%length

       !    write (output_unit_write,*) params%mpi_params%rank,'RR',RR
       !    write (output_unit_write,*) params%mpi_params%rank,'ZZ',spp%vars%Y(:,3)
       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)    
       

       ! Weighting parallel velocity

       !    write (output_unit_write,*) params%mpi_params%rank,'vpll',vpll

       !   do pp=1_idef,spp%ppp
       ! NGP weighting
       rind=FLOOR((rm(cc)-dr/2)/dr)+2_ip
       Vpart(rind)=Vpart(rind)+real(flagCon(cc))*real(flagCol(cc))*vpll(cc)

       ! First-order weighting
!       rind=FLOOR(rm(cc)/dr)+1_ip
!       Vpart(rind)=Vpart(rind)+ &
!            vpll(cc)*(r_1D(rind+1)-rm(cc))/dr
!       Vpart(rind+1)=Vpart(rind+1)+ &
!            vpll(cc)*(rm(cc)-r_1D(rind))/dr

       ! Gaussian weighting
       
!       do ii=1_idef,F%dim_1D
!          arg=MIN((r_1D(ii)-rm(cc))**2._rp/(2._rp*sigr**2._rp),100._rp)
!          Vpart(ii)=Vpart(ii)+1/sqrt(2._rp*C_PI*sigr**2._rp)* &
!               exp(-arg)*vpll(cc)           
!       end do

    end do

    ar=F%AB%a
    ! Calculating density of minor radial annulus    
    do ii=1_idef,F%dim_1D
       ! NGP weighting
       if(ii.eq.1) then
          Vden(ii)=Vpart(ii)/(C_PI*dr**2/4)
       else
          Vden(ii)=Vpart(ii)/(2*C_PI*dr**2*(ii-1))
       end if
       ! First-order weighting
!       if(ii.eq.1) then
!          Vden(ii)=Vpart(ii)/(C_PI*dr**2/3)
!       else
!          Vden(ii)=Vpart(ii)/(2*C_PI*dr**2*(ii-1))
!       end if

       ! Gaussian weighting
!       arg=MIN(r_1D(ii)**2._rp/(2._rp*sigr**2._rp),100._rp)
!       arg1=MIN((ar-r_1D(ii))**2._rp/(2._rp*sigr**2._rp),100._rp)
!       arg2=MIN((ar-r_1D(ii))/(sqrt(2._rp)*sigr),10._rp)
!       arg3=MIN((r_1D(ii))/(sqrt(2._rp)*sigr),10._rp)

!       Ai(ii)=sqrt(C_PI*sigr)*(sqrt(2._rp)*sigr*(exp(-arg)- &
!            exp(-arg1))+r_1D(ii)*sqrt(C_PI)* &
!            (erf(arg2)-erf(-arg3)))
!       Vden(ii)=Vpart(ii)/Ai(ii)

    end do
    
  end subroutine calculate_SC_p