calculate_SC_E1D Subroutine

  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