FIXME : Add documentation
subroutine ginit_recon
use mp, only: proc0
use species, only: spec, has_electron_species
use theta_grid, only: ntgrid, theta
use kt_grids, only: naky, ntheta0, reality
use le_grids, only: forbid
use gs2_save, only: gs2_restore
use dist_fn_arrays, only: g, gnew, vpa, vperp2
use fields_arrays, only: phi, apar, bpar
use fields_arrays, only: phinew, aparnew, bparnew
use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
use file_utils, only: error_unit
use run_parameters, only: has_phi, has_apar, has_bpar
use constants, only: zi
use array_utils, only: copy, zero_array
implicit none
complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phiz, odd
real, dimension (-ntgrid:ntgrid) :: dfac, ufac, tparfac, tperpfac, ct, st, c2t, s2t
integer :: iglo, istatus, ierr
integer :: ig, ik, it, il, is, j
call gs2_restore (g, scale, istatus, has_phi, has_apar, has_bpar)
if (istatus /= 0) then
ierr = error_unit()
if (proc0) write(ierr,*) "Error reading file: ", trim(restart_file)
call zero_array(g)
end if
do iglo = g_lo%llim_proc, g_lo%ulim_proc
ik = ik_idx(g_lo,iglo)
if (ik == 1) cycle
g (:,1,iglo) = 0.
g (:,2,iglo) = 0.
end do
phinew(:,:,2:naky) = 0.
aparnew(:,:,2:naky) = 0.
bparnew(:,:,2:naky) = 0.
phi(:,:,2:naky) = 0.
apar(:,:,2:naky) = 0.
bpar(:,:,2:naky) = 0.
call zero_array(phiz) ; call zero_array(odd)
do j = 1, 2
ik = ikk(j)
it = itt(j)
do ig = -ntgrid, ntgrid
phiz(ig,it,ik) = cmplx(refac, imfac)
end do
end do
odd = zi * phiz
! reality condition for k_theta = 0 component:
if (reality) then
do it = 1, ntheta0/2
phiz(:,it+(ntheta0+1)/2,1) = conjg(phiz(:,(ntheta0+1)/2+1-it,1))
odd (:,it+(ntheta0+1)/2,1) = conjg(odd (:,(ntheta0+1)/2+1-it,1))
enddo
end if
if (even) then
ct = cos(theta)
st = sin(theta)
c2t = cos(2.*theta)
s2t = sin(2.*theta)
else
ct = sin(theta)
st = cos(theta)
c2t = sin(2.*theta)
s2t = cos(2.*theta)
end if
dfac = den0 + den1 * ct + den2 * c2t
ufac = upar0 + upar1* st + upar2* s2t
tparfac = tpar0 + tpar1* ct + tpar2* c2t
tperpfac = tperp0 + tperp1*ct + tperp2*c2t
! charge dependence keeps initial Phi from being too small
do iglo = g_lo%llim_proc, g_lo%ulim_proc
ik = ik_idx(g_lo,iglo)
it = it_idx(g_lo,iglo)
il = il_idx(g_lo,iglo)
is = is_idx(g_lo,iglo)
g(:,1,iglo) = phiinit* &!spec(is)%z* &
( dfac*spec(is)%dens0 * phiz(:,it,ik) &
+ 2.*ufac* vpa(:,1,iglo) * odd (:,it,ik) &
+ tparfac*(vpa(:,1,iglo)**2-0.5) * phiz(:,it,ik) &
+tperpfac*(vperp2(:,iglo)-1.) * phiz(:,it,ik))
where (forbid(:,il)) g(:,1,iglo) = 0.0
g(:,2,iglo) = phiinit* &!spec(is)%z* &
( dfac*spec(is)%dens0 * phiz(:,it,ik) &
+ 2.*ufac* vpa(:,2,iglo) * odd (:,it,ik) &
+ tparfac*(vpa(:,2,iglo)**2-0.5) * phiz(:,it,ik) &
+tperpfac*(vperp2(:,iglo)-1.) * phiz(:,it,ik))
where (forbid(:,il)) g(:,2,iglo) = 0.0
! if (il == ng2+1) g(:,:,iglo) = 0.0
end do
! if (has_electron_species(spec)) then
! call flae (g, gnew)
! g = g - gnew
! end if
call copy(g, gnew)
end subroutine ginit_recon