ginit_recon Subroutine

private subroutine ginit_recon()

FIXME : Add documentation

Arguments

None

Contents

Source Code


Source Code

  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