FIXME : Add documentation
DD>Make sure gamtots are single valued
subroutine init_fieldeq
use dist_fn_arrays, only: aj0, aj1, vperp2, g_work
use species, only: nspec, spec, has_electron_species, has_ion_species, nonmaxw_corr
use theta_grid, only: ntgrid
use kt_grids, only: naky, ntheta0, aky, kperp2
use le_grids, only: integrate_species
use gs2_layouts, only: g_lo, ie_idx, is_idx
use run_parameters, only: tite
implicit none
integer :: iglo, isgn
integer :: ik, it, ie, is
complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: tot
real, dimension (nspec) :: wgt
if (feqinit) return
feqinit = .true.
allocate (gridfac1(-ntgrid:ntgrid,ntheta0,naky))
gridfac1 = 1.0
select case (boundary_option_switch)
case (boundary_option_self_periodic)
! nothing
case (boundary_option_linked)
do it = 1, ntheta0
do ik = 1, naky
if (aky(ik) == 0.0) cycle
if (itleft(it, ik) < 0) gridfac1(-ntgrid,it,ik) = gridfac
if (itright(it, ik) < 0) gridfac1(ntgrid,it,ik) = gridfac
end do
end do
case default
do ik = 1, naky
if (aky(ik) == 0.0) cycle
gridfac1(-ntgrid,:,ik) = gridfac
gridfac1(ntgrid,:,ik) = gridfac
end do
end select
allocate (gamtot(-ntgrid:ntgrid,ntheta0,naky))
allocate (gamtot1(-ntgrid:ntgrid,ntheta0,naky))
allocate (gamtot2(-ntgrid:ntgrid,ntheta0,naky))
if (adiabatic_option_switch == adiabatic_option_fieldlineavg) then
allocate (gamtot3(-ntgrid:ntgrid,ntheta0,naky))
endif
do iglo = g_lo%llim_proc, g_lo%ulim_proc
ie = ie_idx(g_lo,iglo)
is = is_idx(g_lo,iglo)
do isgn = 1, 2
g_work(:,isgn,iglo) = (1.0 - aj0(:,iglo)**2)*nonmaxw_corr(ie,is)
end do
end do
wgt = spec%z*spec%z*spec%dens/spec%temp
call integrate_species (g_work, wgt, tot)
gamtot = real(tot) + kperp2*poisfac
do iglo = g_lo%llim_proc, g_lo%ulim_proc
ie = ie_idx(g_lo,iglo)
is = is_idx(g_lo,iglo)
do isgn = 1, 2
g_work(:,isgn,iglo) = aj0(:,iglo)*aj1(:,iglo) &
*2.0*vperp2(:,iglo) * nonmaxw_corr(ie,is)
end do
end do
wgt = spec%z*spec%dens
call integrate_species (g_work, wgt, tot)
gamtot1 = real(tot)
do iglo = g_lo%llim_proc, g_lo%ulim_proc
ie = ie_idx(g_lo,iglo)
is = is_idx(g_lo,iglo)
do isgn = 1, 2
g_work(:,isgn,iglo) = aj1(:,iglo)**2*2.0*vperp2(:,iglo)**2*nonmaxw_corr(ie,is)
end do
end do
wgt = spec%temp*spec%dens
call integrate_species (g_work, wgt, tot)
gamtot2 = real(tot)
!<DD>Make sure gamtots are single valued
if(esv)then
call ensure_single_val_fields_pass_r(gamtot)
call ensure_single_val_fields_pass_r(gamtot1)
call ensure_single_val_fields_pass_r(gamtot2)
endif
! adiabatic electrons
if (.not. has_electron_species(spec) .or. .not. has_ion_species(spec) ) then
if (adiabatic_option_switch == adiabatic_option_yavg) then
do ik = 1, naky
if (aky(ik) > epsilon(0.0)) gamtot(:,:,ik) = gamtot(:,:,ik) + tite
end do
elseif (adiabatic_option_switch == adiabatic_option_fieldlineavg) then
gamtot = gamtot + tite
gamtot3 = (gamtot-tite) / gamtot
where (gamtot3 < 2.*epsilon(0.0)) gamtot3 = 1.0
else
gamtot = gamtot + tite
endif
endif
end subroutine init_fieldeq