subroutine calculate_SC_E1D(params,F,Vden)
TYPE(FIELDS), INTENT(INOUT) :: F
TYPE(KORC_PARAMS), INTENT(IN) :: params
real(rp),dimension(F%dim_1D),intent(in) :: Vden
real(rp),dimension(F%dim_1D) :: Jsamone,Jsamall,Jexp,dJdt
real(rp),dimension(F%dim_1D) :: a,b,c,u,gam,r
real(rp) :: dr,bet
integer :: ii
INTEGER :: mpierr
! if (params%mpi_params%rank .EQ. 0) then
! write(output_unit_write,*) 'Calculating SC_E1D'
! end if
dr=F%r_1D(2)-F%r_1D(1)
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:5)
Jexp=Jsamall*F%Ip0
F%J3_SC_1D%PHI=F%J2_SC_1D%PHI
F%J2_SC_1D%PHI=F%J1_SC_1D%PHI
F%J1_SC_1D%PHI=Jexp
! Calculating time-derivative of E_phi
dJdt=(3*F%J1_SC_1D%PHI-4*F%J2_SC_1D%PHI+F%J3_SC_1D%PHI)/ &
(2*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*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 calculate_SC_E1D