FIXME : Add documentation
subroutine init_source_term
use dist_fn_arrays, only: vpac, aj0
use run_parameters, only: has_apar
use gs2_time, only: wunits, code_dt
use species, only: spec,nspec,nonmaxw_corr
use hyper, only: D_res
use theta_grid, only: ntgrid, itor_over_b
use le_grids, only: negrid, energy
use constants, only: zi,pi
use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, ie_idx, is_idx
#ifdef LOWFLOW
use lowflow, only: vparterm, wdfac, wstarfac
#endif
implicit none
integer :: iglo
integer :: ig, ik, it, il, ie, is, isgn
!Initialise ufac for use in set_source
if (.not. allocated(ufac)) then
allocate (ufac(negrid, nspec))
do ie = 1, negrid
do is = 1, nspec
ufac(ie, is) = (2.0*spec(is)%uprim &
+ spec(is)%uprim2*energy(ie,is)**(1.5)*sqrt(pi)/4.0)
end do
end do
endif
if(.not.opt_source) return
!Setup the source coefficients
!See comments in get_source_term and set_source
!for information on these terms.
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(iglo, ie, il, ik, it, is, isgn, ig) &
!$OMP SHARED(g_lo, ntgrid, source_coeffs, vpar, &
#ifdef LOWFLOW
!$OMP wdfac, wstarfac, vparterm, &
#endif
!$OMP vpac, code_dt, wunits, ufac, omprimfac, g_exb, itor_over_B, &
!$OMP spec, has_apar, aj0, D_res, nonmaxw_corr, wstar, wdrift) &
!$OMP SCHEDULE(static)
do iglo=g_lo%llim_proc,g_lo%ulim_proc
ie=ie_idx(g_lo,iglo)
il=il_idx(g_lo,iglo)
ik=ik_idx(g_lo,iglo)
it=it_idx(g_lo,iglo)
is=is_idx(g_lo,iglo)
do isgn=1,2
do ig=-ntgrid,ntgrid-1
#ifdef LOWFLOW
!Phi m
source_coeffs(1,ig,isgn,iglo)=-2*vparterm(ig,isgn,iglo)
!Phi p
source_coeffs(2,ig,isgn,iglo)=-zi*wdfac(ig,isgn,iglo) &
+ zi*(wstarfac(ig,isgn,iglo) &
+ vpac(ig,isgn,iglo)*code_dt*wunits(ik)*ufac(ie,is) &
-2.0*omprimfac*vpac(ig,isgn,iglo)*code_dt*wunits(ik)*g_exb*itor_over_B(ig)/spec(is)%stm)
if(has_apar)then
!Apar m
source_coeffs(3,ig,isgn,iglo)=-spec(is)%zstm*&
vpac(ig,isgn,iglo)*(aj0(ig+1,iglo)+aj0(ig,iglo))*0.5
!Apar p
source_coeffs(4,ig,isgn,iglo)=D_res(it,ik)*spec(is)%zstm*&
vpac(ig,isgn,iglo)-spec(is)%stm*vpac(ig,isgn,iglo)*&
zi*(wstarfac(ig,isgn,iglo) &
+ vpac(ig,isgn,iglo)*code_dt*wunits(ik)*ufac(ie,is) &
-2.0*omprimfac*vpac(ig,isgn,iglo)*code_dt*wunits(ik)*g_exb*itor_over_B(ig)/spec(is)%stm)
endif
#else
!Phi m
source_coeffs(1,ig,isgn,iglo)=-2*vpar(ig,isgn,iglo)*nonmaxw_corr(ie,is)
!Phi p
source_coeffs(2,ig,isgn,iglo)=-zi*wdrift(ig,isgn,iglo)*nonmaxw_corr(ie,is) &
+ zi*(wstar(ik,ie,is) &
+ vpac(ig,isgn,iglo)*code_dt*wunits(ik)*ufac(ie,is) &
-2.0*omprimfac*vpac(ig,isgn,iglo)*nonmaxw_corr(ie,is)*code_dt*wunits(ik)*g_exb*itor_over_B(ig)/spec(is)%stm)
if(has_apar)then
!Apar m
source_coeffs(3,ig,isgn,iglo)=-spec(is)%zstm*&
vpac(ig,isgn,iglo)*nonmaxw_corr(ie,is)*(aj0(ig+1,iglo)+aj0(ig,iglo))*0.5
!Apar p
source_coeffs(4,ig,isgn,iglo)=D_res(it,ik)*spec(is)%zstm*vpac(ig,isgn,iglo)*nonmaxw_corr(ie,is)&
-spec(is)%stm*vpac(ig,isgn,iglo)*&
zi*(wstar(ik,ie,is) &
+ vpac(ig,isgn,iglo)*code_dt*wunits(ik)*ufac(ie,is) &
-2.0*omprimfac*vpac(ig,isgn,iglo)*code_dt*wunits(ik)*g_exb*itor_over_B(ig)/spec(is)%stm)
endif
#endif
enddo
enddo
enddo
!$OMP END PARALLEL DO
end subroutine init_source_term