Energy diffusion subroutine used with energy layout (not le_layout) this is always the case when initializing the conserving terms, otherwise is the case if use_le_layout is no specified in the input file.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| complex, | intent(inout), | dimension (-ntgrid:,:,g_lo%llim_proc:) | :: | g |
subroutine solfp_ediffuse_standard_layout (g)
use theta_grid, only: ntgrid
use le_grids, only: negrid, forbid, energy_map
use gs2_layouts, only: ig_idx, it_idx, ik_idx, il_idx, is_idx, e_lo, g_lo
use redistribute, only: gather, scatter
use run_parameters, only: ieqzip
use kt_grids, only: kwork_filter
use array_utils, only: zero_array
implicit none
complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in out) :: g
integer :: ie, is, ig, il, it, ik, ielo
complex, dimension (negrid) :: delta
complex, dimension (:,:), allocatable :: ged
allocate (ged(negrid+1,e_lo%llim_proc:e_lo%ulim_alloc)) ; call zero_array(ged)
call gather (energy_map, g, ged)
! solve for ged row by row
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(ielo, it, ik, is, ig, il, delta, ie) &
!$OMP SHARED(e_lo, kwork_filter, ieqzip, vnew, force_collisions, forbid, negrid, &
!$OMP ged, eql, ec1, ebetaa) &
!$OMP SCHEDULE(static)
do ielo = e_lo%llim_proc, e_lo%ulim_proc
is = is_idx(e_lo, ielo)
ik = ik_idx(e_lo, ielo)
if ((abs(vnew(ik, 1, is)) < 2.0 * epsilon(0.0)) .and. .not. force_collisions) cycle
it = it_idx(e_lo, ielo)
if (kwork_filter(it, ik))cycle
if (ieqzip(it, ik)) cycle
ig = ig_idx(e_lo, ielo)
il = il_idx(e_lo, ielo)
if (forbid(ig, il)) cycle
delta(1) = ged(1, ielo)
do ie = 1, negrid-1
delta(ie+1) = ged(ie+1, ielo) - eql(ie+1, ielo) * delta(ie)
end do
ged(negrid+1, ielo) = 0.0
do ie = negrid, 1, -1
ged(ie, ielo) = (delta(ie) - ec1(ie, ielo) * ged(ie+1, ielo)) * ebetaa(ie, ielo)
end do
end do
!$OMP END PARALLEL DO
call scatter (energy_map, ged, g)
deallocate (ged)
end subroutine solfp_ediffuse_standard_layout