init_SC_E1D Subroutine

  subroutine init_SC_E1D(params,F,spp)
 
    TYPE(FIELDS), INTENT(INOUT)                 :: F
    TYPE(KORC_PARAMS), INTENT(IN) 		:: params
    TYPE(SPECIES), INTENT(IN)    :: spp
    real(rp),dimension(F%dim_1D) :: Vpart
    real(rp),dimension(spp%ppp) :: RR,ZZ,rm,vpll
    real(rp),dimension(F%dim_1D) :: Vden,Jsamone,Jsamall,Jexp,dJdt
    real(rp),dimension(F%dim_1D) :: a,b,c,u,gam,r,r_1D,Ai
    real(rp) :: dr,Isam,bet,sigr,ar,arg,arg1,arg2,arg3
    integer :: pp,ii,rind
    INTEGER 				:: mpierr

!    if (params%mpi_params%rank .EQ. 0) then
!       write(output_unit_write,*) 'Calculating SC_E1D'
!    end if

    ! 1D nearest grid point weighting in minor radius
    
    RR=spp%vars%Y(:,1)
    ZZ=spp%vars%Y(:,3)
    rm=sqrt((RR-F%Ro)**2._rp+(ZZ-F%Zo)**2._rp)*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
    
    
    dr=F%r_1D(2)-F%r_1D(1)

    vpll=spp%vars%V(:,1)/spp%vars%g

    ! Weighting parallel velocity

!    write (output_unit_write,*) 'vpll',vpll
    
    Vpart=0._rp
    r_1D=F%r_1D
    sigr=dr
    
    do pp=1_idef,spp%ppp
       ! NGP weighting
       rind=FLOOR((rm(pp)-dr/2)/dr)+2_ip
       Vpart(rind)=Vpart(rind)+vpll(pp)

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

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

    ! Calculating density of minor radial annulus

    ar=F%AB%a
    
    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
    
    Jsamone=C_E*Vden

    ! Add sampled current densities from all MPI processes Jsamone,
    ! and output of total sampled current density Jsamall to each
    ! MPI process.
    
    call MPI_ALLREDUCE(Jsamone,Jsamall,F%dim_1D,MPI_REAL8,MPI_SUM, &
         MPI_COMM_WORLD,mpierr)
    
!    write(output_unit_write,*) 'Jsam: ',Jsamall(1:10)

    ! Integrating current density to scale total current to
    ! experimentally determined total current


    Isam=0._rp
    do ii=1_idef,F%dim_1D
       if ((ii.eq.1).or.(ii.eq.F%dim_1D)) then
          Isam=Isam+Jsamall(ii)*r_1D(ii)/2._rp
       else 
          Isam=Isam+Jsamall(ii)*r_1D(ii)
       end if
    end do
    Isam=2._rp*C_PI*Isam*dr
!    write(output_unit_write,*) params%mpi_params%rank,'Isam: ',Isam

    F%Ip0=F%Ip_exp/Isam

    
    Jexp=Jsamall*F%Ip0

    F%J3_SC_1D%PHI=Jexp 
    F%J2_SC_1D%PHI=Jexp
    F%J1_SC_1D%PHI=Jexp        
    
    ! Calculating time-derivative of E_phi
    
    dJdt=(3._rp*F%J1_SC_1D%PHI-4._rp*F%J2_SC_1D%PHI+F%J3_SC_1D%PHI)/ &
         (2._rp*F%dt_E_SC)
    
!    write(output_unit_write,*) params%mpi_params%rank,'J(1)',F%J_SC_1D%PHI(1)
    
    ! Solving 1D Poisson equation with tridiagonal matrix solve

    a=0._rp
    b=-2._rp
    c=0._rp
    u=0._rp
    gam=0._rp
!    r=-2*dr**2*C_MU*Jexp
    r=2*dr**2*C_MU*dJdt

    do ii=2_idef,F%dim_1D
       a(ii)=(REAL(ii)-2._rp)/(REAL(ii)-1._rp)
       c(ii)=REAL(ii)/(REAL(ii)-1._rp)
    end do

    bet=b(2)
    u(2)=r(2)/bet
    do ii=3_idef,F%dim_1D-1
       gam(ii)=c(ii-1)/bet
       bet=b(ii)-a(ii)*gam(ii)
       if (bet.eq.0) then
          stop 'tridiag failed'         
       end if
       u(ii)=(r(ii)-a(ii)*u(ii-1))/bet
    end do

    do ii=F%dim_1D-2,2,-1
       u(ii)=u(ii)-gam(ii+1)*u(ii+1)
    end do

    u(1)=(4._rp*u(2)-u(3))/3._rp

    ! Writing over F%A* data
    
!    F%A3_SC_1D%PHI=F%A2_SC_1D%PHI
!    F%A2_SC_1D%PHI=F%A1_SC_1D%PHI
!    F%A1_SC_1D%PHI=u
    
!    if (init) then
!       F%A3_SC_1D%PHI=F%A1_SC_1D%PHI
!       F%A2_SC_1D%PHI=F%A1_SC_1D%PHI 
!    end if

!    write(output_unit_write,*) params%mpi_params%rank,'A1(1)',F%A1_SC_1D%PHI(1)
!    write(output_unit_write,*) params%mpi_params%rank,'A2(1)',F%A2_SC_1D%PHI(1)
!    write(output_unit_write,*) params%mpi_params%rank,'A3(1)',F%A3_SC_1D%PHI(1)
    
    ! Calculating inductive E_phi

!    F%E_SC_1D%PHI=-(3*F%A1_SC_1D%PHI-4*F%A2_SC_1D%PHI+F%A3_SC_1D%PHI)/ &
!         (2*F%dt_E_SC)

    F%E_SC_1D%PHI=u

    if (params%mpi_params%rank.eq.0) then
       write(output_unit_write,*) 'J1(2)',F%J1_SC_1D%PHI(2)
       write(output_unit_write,*) 'J2(2)',F%J2_SC_1D%PHI(2)
       write(output_unit_write,*) 'J3(2)',F%J3_SC_1D%PHI(2)
       
       write(output_unit_write,*) 'E(1)',F%E_SC_1D%PHI(1)
    end if
       
    ! Normalizing inductive E_phi
    
    F%E_SC_1D%PHI=F%E_SC_1D%PHI/params%cpp%Eo

#ifdef PSPLINE
    call initialize_SC1D_field_interpolant(params,F)
#endif
    
  end subroutine init_SC_E1D