subroutine init_SC_E1D_FS(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) :: dintJphidPSIP,PSIP_1D
real(rp),dimension(spp%ppp) :: PSIP,vpll
real(rp),dimension(F%dim_1D) :: Jsamall,Jexp,dJdt
real(rp),dimension(F%dim_1D) :: a,b,c,u,gam,r,alpha,beta,gamma
real(rp) :: dPSIP,Isam,bet,arg,arg1,PSIP_lim,sigPSIP
integer :: pp,ii,PSIPind
INTEGER :: mpierr
PSIP_1D=F%PSIP_1D
dPSIP=PSIP_1D(2)-PSIP_1D(1)
PSIP=spp%vars%PSI_P*(params%cpp%Bo*params%cpp%length**2)
sigPSIP=dPSIP
vpll=spp%vars%V(:,1)/spp%vars%g
dintJphidPSIP=0._rp
do pp=1_idef,spp%ppp
if (PSIP(pp).lt.0._rp) PSIP(pp)=0._rp
PSIPind=FLOOR(PSIP(pp)/dPSIP)+1_ip
! NGP weighting
! dintJphidPSIP(PSIPind)=dintJphidPSIP(PSIPind)+vpll(pp)
! First-order weighting
! dintJphidPSIP(PSIPind)=dintJphidPSIP(PSIPind)+ &
! vpll(pp)*(PSIP_1D(PSIPind+1)-PSIP(pp))/dPSIP
! dintJphidPSIP(PSIPind+1)=dintJphidPSIP(PSIPind+1)+ &
! vpll(pp)*(PSIP(pp)-PSIP_1D(PSIPind))/dPSIP
! write(output_unit_write,*) PSIP(pp),PSIP_1D(PSIPind),dPSIP
! Gaussian weighting
do ii=1_idef,F%dim_1D
arg=MIN((PSIP_1D(ii)-PSIP(pp))**2._rp/ &
(2._rp*sigPSIP**2._rp),100._rp)
dintJphidPSIP(ii)=dintJphidPSIP(ii)+ &
vpll(pp)*exp(-arg)
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
! Add sampled current densities from all MPI processes Jsamone,
! and output of total sampled current density Jsamall to each
! MPI process.
call MPI_ALLREDUCE(dintJphidPSIP,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)/2._rp
else
Isam=Isam+Jsamall(ii)
end if
end do
Isam=Isam*dPSIP
! 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
alpha=F%ddMagPsiSqdPsiPSq
beta=F%dMagPsiSqdPsiP
gamma=C_MU*dJdt
a=-alpha*dPSIP/2._rp+beta
b=-2._rp*beta
c=alpha*dPSIP/2._rp+beta
u=0._rp
gam=0._rp
! r=-2*dr**2*C_MU*Jexp
r=dPSIP**2*gamma
c(2)=c(2)-a(2)*a(1)/c(1)
b(2)=b(2)-a(2)*b(1)/c(1)
r(2)=r(2)-a(2)*r(1)/c(1)
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)=2*u(2)-u(3)
F%E_SC_1D%PHI=u
if (params%mpi_params%rank.eq.0) then
write(output_unit_write,*) 'J1(1)',F%J1_SC_1D%PHI(1)
write(output_unit_write,*) 'J2(1)',F%J2_SC_1D%PHI(1)
write(output_unit_write,*) 'J3(1)',F%J3_SC_1D%PHI(1)
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_FS(params,F)
#endif
end subroutine init_SC_E1D_FS