1970319704

Committed 06 Aug 2025 03:50PM UTC coverage: 8.18% (+0.002%) from 8.178%

Build # 1970319704

Build Type

push

gitlab-ci

Committed by David Dickinson

Commit Message

Merged in bugfix/fix_issue_88_multiple_ncheck (pull request #1109)

Run Details

3667 of 44831 relevant lines covered (8.18%)

124503.35 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/src/diagnostics/diagnostics_velocity_space.fpp

!> A module which writes out quantities which writes out 
!! diagnostic quantities which assess whether the velocity 
!! space resolution is sufficient.
module diagnostics_velocity_space
  implicit none
  private

  public :: init_diagnostics_velocity_space, write_velocity_space_checks, get_gtran
  public :: write_collision_error, finish_diagnostics_velocity_space, get_verr
  public :: collision_error, finish_legendre_transform, setup_legendre_transform
  public :: init_weights, finish_weights

  ! GTRAN
  real, dimension (:,:), allocatable :: lpl
  real, dimension (:,:,:), allocatable :: lpe
  real, dimension (:,:,:,:), allocatable :: lpt
  ! VERR
  real, dimension (:,:,:), allocatable :: wlmod, wtmod, wmod

contains

  !> FIXME : Add documentation
  subroutine init_diagnostics_velocity_space(gnostics)
    use diagnostics_config, only: diagnostics_type
    use collisions, only: collision_model_switch, collision_model_lorentz
    use collisions, only: init_lorentz_error, collision_model_lorentz_test
    implicit none
    type(diagnostics_type), intent(inout) :: gnostics 
    
    if (.not. gnostics%write_verr) return
    
    ! initialize weights for less accurate integrals used
    ! to provide an error estimate for v-space integrals (energy and untrapped)
    call init_weights
    call setup_legendre_transform
    
    if (gnostics%write_cerr) then
       if (collision_model_switch == collision_model_lorentz .or. &
            collision_model_switch == collision_model_lorentz_test) then
          call init_lorentz_error
       else
          gnostics%write_cerr = .false.
       end if
    end if
  end subroutine init_diagnostics_velocity_space

  !> FIXME : Add documentation  
  subroutine finish_diagnostics_velocity_space()
    implicit none
    call finish_weights
    call finish_legendre_transform
  end subroutine finish_diagnostics_velocity_space

  !> FIXME : Add documentation  
  subroutine write_velocity_space_checks(gnostics, vary_vnew_only)
    use mp, only: proc0
    use le_grids, only: grid_has_trapped_particles
    use fields_arrays, only: phinew, bparnew
    use gs2_time, only: user_time
    use collisions, only: vnmult
    use species, only: spec
    use diagnostics_config, only: diagnostics_type
#ifdef NETCDF
    use gs2_io, only: starts, time_dim, dim_2, dim_3, dim_5
    use neasyf, only: neasyf_write
#endif
    implicit none
    type(diagnostics_type), intent(in) :: gnostics
    logical, intent(in) :: vary_vnew_only
    real, dimension (5,2) :: errest
    integer, dimension (5,3) :: erridx
    real :: geavg, glavg, gtavg

    errest = 0.0; erridx = 0
    geavg = 0.0 ; glavg = 0.0 ; gtavg = 0.0

    ! error estimate obtained by comparing standard integral with less-accurate integral
    call get_verr (errest, erridx, phinew, bparnew)

    if (vary_vnew_only) return

    ! error estimate based on monitoring amplitudes of legendre polynomial coefficients
    call get_gtran (geavg, glavg, gtavg, phinew, bparnew)

#ifdef NETCDF
    if (gnostics%writing) then
       call neasyf_write(gnostics%file_id, "vspace_lpcfrac", [geavg, glavg, gtavg], &
            dim_names=[dim_3, time_dim], start=starts(2, gnostics%nout), &
            long_name="Fraction of free energy contained in the high order coefficients of &
            & the Legendre polynomial transform of (1) energy space, (2) untrapped &
            & pitch angles and (3) trapped pitch angles (each should ideally be < 0.1).  &
            & Note that there are no trapped pitch angles for certain geometries")
       call neasyf_write(gnostics%file_id, "vspace_err", errest, &
            dim_names=[dim_5, dim_2, time_dim], start=starts(3, gnostics%nout), &
            long_name="Estimate of the (1) absolute and (2) relative errors resulting from &
            & velocity space integrals in the calculation of the following quantities &
            & in the given dimensions: (1) k phi, energy (2) k phi, untrapped pitch angles &
            & (3) k phi, trapped pitch angles, (4) k apar, energy, (5) k apar, untrapped &
            & angles. Relative errors should be < 0.1. ")
       call neasyf_write(gnostics%file_id, "vspace_vnewk", &
            [vnmult(1)*spec(1)%vnewk, vnmult(2)*spec(1)%vnewk], &
            dim_names=[dim_2, time_dim], start=starts(2, gnostics%nout), &
            long_name="If the simulation is set to vary the collisionality in order to keep &
            & error in velocity integrals to acceptable levels, contains species 1 &
            & collisionality in (1) pitch angle and (2) energy  ")

       if (gnostics%write_max_verr) then
          call neasyf_write(gnostics%file_id, "vspace_err_maxindex", erridx, &
               dim_names=[dim_5, dim_3, time_dim], start=starts(3, gnostics%nout), &
               long_name="Gives the (1) theta index, (2) ky index and (3) kx index of the maximum &
               & error resulting from the &
               & velocity space integrals in the calculation of the following quantities &
               & in the given dimensions: (1) k phi, energy (2) k phi, untrapped pitch angles &
               & (3) k phi, trapped pitch angles, (4) k apar, energy, (5) k apar, untrapped &
               & angles. Relative errors should be < 0.1. ")
       end if
    end if
#endif
    if (proc0 .and. gnostics%ascii_files%write_to_vres .and. (.not. gnostics%create)) call write_ascii
  contains
    !> FIXME : Add documentation
    subroutine write_ascii
      if (grid_has_trapped_particles()) then
         write(gnostics%ascii_files%lpc,"(4(1x,e13.6))") user_time, geavg, glavg, gtavg
      else
         write(gnostics%ascii_files%lpc,"(3(1x,e13.6))") user_time, geavg, glavg
      end if
      write(gnostics%ascii_files%vres,"(8(1x,e13.6))") user_time, errest(1,2), errest(2,2), errest(3,2), &
           errest(4,2), errest(5,2), vnmult(1)*spec(1)%vnewk, vnmult(2)*spec(1)%vnewk
      if (gnostics%write_max_verr) then
         write(gnostics%ascii_files%vres2,"(3(i8),(1x,e13.6),3(i8),(1x,e13.6),3(i8),(1x,e13.6),3(i8),(1x,e13.6),3(i8),(1x,e13.6))") &
              erridx(1,1), erridx(1,2), erridx(1,3), errest(1,1), &
              erridx(2,1), erridx(2,2), erridx(2,3), errest(2,1), &
              erridx(3,1), erridx(3,2), erridx(3,3), errest(3,1), &
              erridx(4,1), erridx(4,2), erridx(4,3), errest(4,1), &
              erridx(5,1), erridx(5,2), erridx(5,3), errest(5,1)
      end if
    end subroutine write_ascii
  end subroutine write_velocity_space_checks

  !> Calculate and write the collision error
  subroutine write_collision_error (gnostics)
    use diagnostics_config, only: diagnostics_type
    use fields_arrays, only: phinew, bparnew
    type(diagnostics_type), intent(in) :: gnostics
    if (gnostics%write_ascii) call collision_error(gnostics%ascii_files%cres, phinew, bparnew)
  end subroutine write_collision_error

  !> Calculate the collision error and write to text file `<runname>.cres`
  subroutine collision_error(unit, phi, bpar)
    use mp, only: proc0, iproc, send, receive
    use le_grids, only: ng2, jend, nlambda, lambda_map
    use theta_grid, only: ntgrid
    use dist_fn_arrays, only: gnew, g_adjust, g_work, to_g_gs2, from_g_gs2
    use gs2_layouts, only: lz_lo, ig_idx, idx_local, proc_id
    use gs2_layouts, only: ik_idx, ie_idx, is_idx, it_idx, il_idx
    use collisions, only: dtot, fdf, fdb
    use redistribute, only: gather
    use gs2_time, only: user_time
    use array_utils, only: zero_array
    implicit none
    integer, intent(in) :: unit
    complex, dimension(-ntgrid:, :, :), intent(in) :: phi, bpar
    integer :: je, te, ig, il, ip, ilz, ie, is, ik, it
    integer :: igmax, ikmax, itmax, iemax, ilmax, ismax, proc_idx
    complex, dimension (:), allocatable :: ltmp, ftmp
    complex, dimension (:,:), allocatable :: lcoll, fdcoll, glze
    real :: etmp, emax, etot, eavg, edenom, ltmax

    allocate (ltmp(2*nlambda), ftmp(2*nlambda))
    allocate (lcoll(2*nlambda,lz_lo%llim_proc:lz_lo%ulim_alloc))
    allocate (fdcoll(2*nlambda,lz_lo%llim_proc:lz_lo%ulim_alloc))
    allocate (glze(max(2*nlambda,2*ng2+1),lz_lo%llim_proc:lz_lo%ulim_alloc))

    call zero_array(lcoll); call zero_array(fdcoll); call zero_array(glze)
    ltmp = 0.0; ftmp = 0.0
    etmp = 0.0; emax = 0.0; etot = 0.0; eavg = 0.0; edenom = 1.0; ltmax = 1.0

    ! convert gnew from g to h
    call g_adjust(gnew, phi, bpar, direction = from_g_gs2)
    g_work = gnew
    ! convert gnew from h back to g
    call g_adjust(gnew, phi, bpar, direction = to_g_gs2)

    ! map from g_work(ig,isgn,iglo) to glze(il,ilz)
    ! Note lambda_map is only defined if use_le_layout = F (not the default).
    ! We force write_cerr = F if use_le_layout = T
    call gather (lambda_map, g_work, glze)

    ! loop over ig, isign, ik, it, ie, is
    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ilz, ig, je, te, il, ip) &
    !$OMP SHARED(lz_lo, jend, ng2, dtot, glze, fdf, fdb, fdcoll) &
    !$OMP REDUCTION(+: lcoll) &
    !$OMP SCHEDULE(static)
    do ilz = lz_lo%llim_proc, lz_lo%ulim_proc
       ig = ig_idx(lz_lo,ilz)

       if (jend(ig) == 0) then      ! no trapped particles
          ! je = number of + vpa grid pts
          ! te = number of + and - vpa grid pts
          je = ng2
          te = 2*ng2

          ! find d/d(xi) ((1+xi**2)( d g(xi)/ d(xi) )) at each xi
          ! using lagrange (lcoll) and finite difference (fdcoll)
          il = 1
          do ip = il, il+2
             lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,il,ip)*glze(ip,ilz)
          end do

          il = 2
          do ip = il-1, il+1
             lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,il,ip-il+2)*glze(ip,ilz)
          end do

          do il=3,ng2
             do ip=il-2,il+2
                lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,il,ip-il+3)*glze(ip,ilz)
             end do
          end do

          do il=ng2+1, 2*ng2-2
             do ip = il-2,il+2
                lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,2*ng2-il+1,il-ip+3)*glze(ip,ilz)
             end do
          end do

          il = 2*ng2-1
          do ip = il-1, il+1
             lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,2,il-ip+2)*glze(ip,ilz)
          end do

          il = 2*ng2
          do ip = il-2, il
             lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,1,il-ip+1)*glze(ip,ilz)
          end do

          ! deal with xi from 1-eps -> eps
          do il=2,ng2
             fdcoll(il,ilz) = fdf(ig,il)*(glze(il+1,ilz) - glze(il,ilz)) - fdb(ig,il)*(glze(il,ilz) - glze(il-1,ilz))
          end do
          ! deal with xi from -eps -> -1+eps
          do il=ng2+1, 2*ng2-1
             fdcoll(il,ilz) = fdb(ig,2*ng2-il+1)*(glze(il+1,ilz) - glze(il,ilz)) - fdf(ig,2*ng2-il+1)*(glze(il,ilz) - glze(il-1,ilz))
          end do

          fdcoll(1,ilz) = fdf(ig,1)*(glze(2,ilz) - glze(1,ilz))
          fdcoll(2*ng2,ilz) = -fdf(ig,1)*(glze(2*ng2,ilz) - glze(2*ng2-1,ilz))

       else       ! trapped particle runs          
          je = jend(ig)
          te = 2*je - 1

          il = 1
          do ip = il, il+2
             lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,il,ip)*glze(ip,ilz)
          end do

          il = 2
          do ip = il-1, il+1
             lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,il,ip-il+2)*glze(ip,ilz)
          end do

          do il=3,je
             do ip=il-2,il+2
                lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,il,ip-il+3)*glze(ip,ilz)
             end do
          end do

          do il=je+1, te-2
             do ip = il-2,il+2
                lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,te-il+1,il-ip+3)*glze(ip,ilz)
             end do
          end do

          il = te-1
          do ip = il-1, il+1
             lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,2,il-ip+2)*glze(ip,ilz)
          end do

          il = te
          do ip = il-2, il
             lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,1,il-ip+1)*glze(ip,ilz)
          end do

          ! is il=je handled correctly here?
          do il=2,je
             fdcoll(il,ilz) = fdf(ig,il)*(glze(il+1,ilz) - glze(il,ilz)) - fdb(ig,il)*(glze(il,ilz) - glze(il-1,ilz))
          end do
          do il=je+1,te-1
             fdcoll(il,ilz) = fdb(ig,te-il+1)*(glze(il+1,ilz) - glze(il,ilz)) - fdf(ig,te-il+1)*(glze(il,ilz) - glze(il-1,ilz))
          end do

          fdcoll(1,ilz) = fdf(ig,1)*(glze(2,ilz) - glze(1,ilz))
          fdcoll(te,ilz) = -fdf(ig,1)*(glze(te,ilz) - glze(te-1,ilz))

       end if
    end do
    !$OMP END PARALLEL DO

    do ilz=lz_lo%llim_world, lz_lo%ulim_world
       ig = ig_idx(lz_lo, ilz)
       ik = ik_idx(lz_lo, ilz)
       it = it_idx(lz_lo, ilz)
       ie = ie_idx(lz_lo, ilz)
       is = is_idx(lz_lo, ilz)
       je = jend(ig)

       if (je == 0) then
          te = 2 * ng2
       else
          te = 2 * je - 1
       end if

       proc_idx = proc_id(lz_lo, ilz)
       if (proc_idx == iproc) then
          if (proc0) then 
             ltmp = lcoll(:,ilz)
             ftmp = fdcoll(:,ilz)
          else
             call send (lcoll(:,ilz), 0)
             call send (fdcoll(:,ilz), 0)
          endif
       else if (proc0) then
          call receive (ltmp, proc_idx)
          call receive (ftmp, proc_idx)
       endif

       ! Only proc0 has ltmp and ftmp set so on that proc calculate the error
       ! (difference) and look for maximum
       if (proc0) then
          do il = 1, te
             etmp = abs(ltmp(il) - ftmp(il))

             if (etmp > emax) then
                emax = etmp
                ltmax = abs(ltmp(il))
                ikmax = ik
                itmax = it
                iemax = ie
                ismax = is
                ilmax = il
                igmax = ig
             end if

             etot = etot + etmp
             edenom = edenom + abs(ltmp(il))
          end do
       end if
    end do

    if (proc0) then
       eavg = etot / edenom
       emax = emax / ltmax
       write(unit,"((1x,e13.6),6(i8),2(1x,e13.6))") user_time, &
            igmax, ikmax, itmax, iemax, ilmax, ismax, emax, eavg
    end if

    deallocate (lcoll, fdcoll, glze, ltmp, ftmp)

  end subroutine collision_error

  !> Error estimate obtained by comparing standard integral with less-accurate integral
  !>
  !> Estimate of the (1) absolute and (2) relative errors resulting from velocity space
  !> integrals in the calculation of the following quantities in the given dimensions: (1)
  !> k phi, energy (2) k phi, untrapped pitch angles (3) k phi, trapped pitch angles, (4)
  !> k apar, energy, (5) k apar, untrapped angles. Relative errors should be < 0.1.
  !>
  !> @note Should this be extended to include error estimates on bpar as well?
  subroutine get_verr (errest, erridx, phi, bpar)
    use le_grids, only: integrate_species, nesub, new_trap_int, grid_has_trapped_particles
    use le_grids, only: ng2, nlambda
    use theta_grid, only: ntgrid
    use kt_grids, only: ntheta0, naky, aky, akx
    use species, only: nspec, spec
    use dist_fn_arrays, only: gnew, aj0, vpa, g_adjust, g_work, to_g_gs2, from_g_gs2
    use run_parameters, only: has_phi, has_apar, beta
    use gs2_layouts, only: g_lo
    use collisions, only: adjust_vnmult
    use array_utils, only: zero_array
    implicit none
    !> Indices of maximum error.
    integer, dimension (5,3), intent (out) :: erridx
    !> Estimated error.
    real, dimension (5,2), intent (out) :: errest
    !> Electrostatic potential and parallel magnetic field
    complex, dimension (-ntgrid:,:,:), intent (in) :: phi, bpar
    integer :: ig, it, ik, iglo, isgn, ntrap
    complex, dimension (:,:,:), allocatable :: phi_app, apar_app
    complex, dimension (:,:,:,:), allocatable :: phi_e, phi_l, phi_t, apar_e, apar_l, apar_t
    real, dimension (:,:), allocatable :: kmax
    real, dimension (:), allocatable :: wgt
    real, dimension(2) :: errtmp
    integer, dimension(3) :: idxtmp
    logical, parameter :: trap_flag = .true.
    logical :: compute_trapped_error

    allocate(wgt(nspec))

    if (has_phi) then
       allocate(phi_app(-ntgrid:ntgrid,ntheta0,naky))
       allocate(phi_e(-ntgrid:ntgrid,ntheta0,naky,nesub))
       allocate(phi_l(-ntgrid:ntgrid,ntheta0,naky,ng2))
    end if

    if (has_apar) then
       allocate(apar_app(-ntgrid:ntgrid,ntheta0,naky))
       allocate(apar_e(-ntgrid:ntgrid,ntheta0,naky,nesub))
       allocate(apar_l(-ntgrid:ntgrid,ntheta0,naky,ng2))
    end if

    compute_trapped_error = grid_has_trapped_particles() .and. new_trap_int

    ! first call to g_adjust converts gyro-averaged dist. fn. (g)
    ! into nonadiabatic part of dist. fn. (h)
    call g_adjust (gnew, phi, bpar, direction = from_g_gs2)

    ! These next if statements and loops are essentially repeating some of getan
    ! and adding in error estimate calls, which are themselves essentially integrate_species
    ! but done multiple times for different v-space weights.

    ! take gyro-average of h at fixed total position (not g.c. position)
    ! Note here phi_app and apar_app are effectively antot and antota that
    ! would be returned from a call to getan_from_dnf(gnew, antot, antota,...)
    if (has_phi) then
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(iglo, isgn, ig) &
       !$OMP SHARED(g_lo, ntgrid, g_work, aj0, gnew) &
       !$OMP COLLAPSE(3) &
       !$OMP SCHEDULE(static)
       do iglo = g_lo%llim_proc, g_lo%ulim_proc
          do isgn = 1, 2
             do ig = -ntgrid, ntgrid
                g_work(ig, isgn, iglo) = aj0(ig, iglo) * gnew(ig, isgn, iglo)
             end do
          end do
       end do
       !$OMP END PARALLEL DO

       wgt = spec%z*spec%dens
       call integrate_species (g_work, wgt, phi_app)

       ! integrates dist fn of each species over v-space
       ! after dropping an energy grid point and returns
       ! phi_e, which contains the integral approximations
       ! to phi for each point dropped
       call eint_error (g_work, wgt, phi_e)

       ! integrates dist fn of each species over v-space
       ! after dropping an untrapped lambda grid point and returns phi_l.
       ! phi_l contains ng2 approximations for the integral over lambda that
       ! come from dropping different pts from the gaussian quadrature grid
       call lint_error (g_work, wgt, phi_l)

       ! next loop gets error estimate for trapped particles, if there are any
       if (compute_trapped_error) then
          ntrap = nlambda - ng2
          allocate(phi_t(-ntgrid:ntgrid, ntheta0, naky, ntrap))
          call zero_array(phi_t)
          call trap_error (g_work, wgt, phi_t)
       end if

    end if

    if (has_apar) then
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(iglo, isgn, ig) &
       !$OMP SHARED(g_lo, ntgrid, g_work, aj0, vpa, gnew) &
       !$OMP COLLAPSE(3) &
       !$OMP SCHEDULE(static)
       do iglo = g_lo%llim_proc, g_lo%ulim_proc
          do isgn = 1, 2
             do ig = -ntgrid, ntgrid
                g_work(ig, isgn, iglo) = aj0(ig, iglo) * vpa(ig, isgn, iglo) * &
                     gnew(ig, isgn, iglo)
             end do
          end do
       end do
       !$OMP END PARALLEL DO

       wgt = 2.0 *beta * spec%z * spec%dens * sqrt(spec%temp / spec%mass)
       call integrate_species (g_work, wgt, apar_app)

       call eint_error (g_work, wgt, apar_e)
       call lint_error (g_work, wgt, apar_l)
       if (compute_trapped_error) then
          ntrap = nlambda - ng2
          allocate(apar_t(-ntgrid:ntgrid, ntheta0, naky, ntrap))
          call zero_array(apar_t)
          call trap_error (g_work, wgt, apar_t)
       end if
    end if
    deallocate (wgt)

    ! second call to g_adjust converts from h back to g
    call g_adjust (gnew, phi, bpar, direction = to_g_gs2)

    allocate (kmax(ntheta0, naky))
    do ik = 1, naky
       do it = 1, ntheta0
          kmax(it,ik) = max(akx(it),aky(ik))
       end do
    end do

    errest = 0.0
    erridx = 0

    if (has_phi) then

       call estimate_error (phi_app, phi_e, kmax, errtmp, idxtmp)
       errest(1,:) = errtmp
       erridx(1,:) = idxtmp

       call estimate_error (phi_app, phi_l, kmax, errtmp, idxtmp)
       errest(2,:) = errtmp
       erridx(2,:) = idxtmp

       if (compute_trapped_error) then
          call estimate_error (phi_app, phi_t, kmax, errtmp, idxtmp, trap_flag)
          errest(3,:) = errtmp
          erridx(3,:) = idxtmp
          deallocate (phi_t)
       end if
       deallocate(phi_app, phi_e, phi_l)
    end if

    ! No trapped errors for apar?
    if (has_apar) then

       call estimate_error (apar_app, apar_e, kmax, errtmp, idxtmp)
       errest(4,:) = errtmp
       erridx(4,:) = idxtmp

       call estimate_error (apar_app, apar_l, kmax, errtmp, idxtmp)
       errest(5,:) = errtmp
       erridx(5,:) = idxtmp

       deallocate(apar_app, apar_e, apar_l)
    end if

    call adjust_vnmult(errest, compute_trapped_error)
  end subroutine get_verr

  !> FIXME : Add documentation
  subroutine estimate_error (app1, app2, kmax_local, errest, erridx, trap)
    use kt_grids, only: naky, ntheta0
    use theta_grid, only: ntgrid
    use le_grids, only: ng2, jend, il_is_trapped
    use optionals, only: get_option_with_default
    implicit none
    complex, dimension (-ntgrid:,:,:), intent (in) :: app1
    complex, dimension(-ntgrid:,:,:,:), intent (in) :: app2
    real, dimension(:, :), intent (in) :: kmax_local
    logical, intent (in), optional :: trap
    real, dimension(2), intent (out) :: errest
    integer, dimension(3), intent (out) :: erridx
    real, dimension(-ntgrid:ntgrid, ntheta0, naky) :: gpavg_arr, gs_arr
    integer, dimension(3) :: max_locs
    integer :: ik, it, ig, end_idx, igmax, ikmax, itmax
    real :: gdsum, gdmax, gnsum
    logical :: do_trap
    do_trap = get_option_with_default(trap, .false.)
    gpavg_arr = 0.0 ; gs_arr = 0.0

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ik, it, ig, end_idx) &
    !$OMP SHARED(naky, ntheta0, ntgrid, do_trap, jend, ng2, app2, &
    !$OMP app1, kmax_local, gpavg_arr, gs_arr) &
    !$OMP COLLAPSE(3) &
    !$OMP SCHEDULE(static)
    do ik = 1, naky
       do it = 1, ntheta0
          do ig = -ntgrid, ntgrid
             if (.not. do_trap .or. il_is_trapped(jend(ig))) then
                if (do_trap) then
                   end_idx = jend(ig) - ng2
                else
                   end_idx = size(app2, 4)
                end if

                gpavg_arr(ig, it, ik) = sum(abs(app1(ig, it, ik) - app2(ig, it, ik, 1:end_idx))) * kmax_local(it, ik) / end_idx
                gs_arr(ig, it, ik) = kmax_local(it, ik) * abs(app1(ig, it, ik))
             end if
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    max_locs = maxloc(gpavg_arr)
    igmax = max_locs(1) - (ntgrid + 1) !Adjust for non-1 lbound
    itmax = max_locs(2)
    ikmax = max_locs(3)

    if (abs(gs_arr(igmax, itmax, ikmax)) < epsilon(0.0)) then
       gdmax = 0.0
    else
       gdmax = gpavg_arr(igmax, itmax, ikmax) / gs_arr(igmax, itmax, ikmax)
    end if

    erridx(1) = igmax
    erridx(2) = ikmax
    erridx(3) = itmax
    errest(1) = gdmax

    gnsum = sum(gpavg_arr)
    gdsum = sum(gs_arr)

    ! Guard against divide by zero. This can commonly occur on the first time step
    ! in runs with fapar /= 0 due to our initial conditions typically ensuring that
    ! the initial apar is identically zero.
    if (abs(gdsum) < epsilon(0.0)) then
       errest(2) = 0.0
    else
       errest(2) = gnsum / gdsum
    end if
  end subroutine estimate_error

  !> Calculates energy grid and (untrapped) pitch angle weights for
  !> integrals which drop a point from the energy/lambda
  !> dimensions. These are only used for estimating the velocity space
  !> error in [[dist_fn::get_verr]] as a part of the adaptive
  !> collisionality algorithm.
  subroutine init_weights
    use file_utils, only: open_output_file, close_output_file
    use constants, only: twopi => dtwopi
    use species, only: nspec, f0_values
    use le_grids, only: get_weights, speed, nesub, ng2, vcut, w, energy, xx, negrid, nmax
    use le_grids, only: grid_has_trapped_particles, wl, nlambda, al, new_trap_int
    use le_grids, only: setup_trapped_lambda_grids_new_trap_int
    use theta_grid, only: bmag, bmax, ntgrid
    implicit none
    real, dimension (:), allocatable :: modzeroes, werrtmp  ! (negrid-2)
    real, dimension (:), allocatable :: lmodzeroes, wlerrtmp ! (ng2-1)
    real, dimension(:,:), allocatable :: dummy, wlerr
    real, dimension(:,:,:), allocatable :: wlterr, werr
    integer :: ipt, ndiv, divmax, is, ntrap, il
    logical :: eflag

    allocate(lmodzeroes(ng2-1), wlerrtmp(ng2-1))
    allocate(wlerr(ng2,ng2))
    wlerr = 0.0; lmodzeroes = 0.0; wlerrtmp = 0.0
    allocate(modzeroes(nesub-1), werrtmp(nesub-1))
    allocate(werr(negrid-1,nesub,nspec))
    werr = 0.0 ; modzeroes = 0.0 ; werrtmp = 0.0

    ! loop to obtain weights for energy grid points.  negrid-1 sets
    ! of weights are needed because we want to compute integrals
    ! for negrid-1 sets of energy points (corresponding to negrid-1
    ! points that we can choose to drop from the guassian quadrature)
    do ipt = 1, nesub
       do is = 1, nspec

         ! drops the point corresponding to ipt from the energy grid
         if (ipt /= 1) modzeroes(:ipt-1) = speed(:ipt-1,is)
         if (ipt /= nesub) modzeroes(ipt:nesub-1) = speed(ipt+1:nesub,is)

         ! get weights for energy grid points
         call get_weights (nmax,0.0,vcut,modzeroes,werrtmp,ndiv,divmax,eflag)

         ! a zero is left in the position corresponding to the dropped point
         if (ipt /= 1) werr(:ipt-1,ipt,is) = werrtmp(:ipt-1)
         if (ipt /= nesub) werr(ipt+1:nesub,ipt,is) = werrtmp(ipt:nesub-1)
         werr(nesub+1:,ipt,is) = w(nesub+1:negrid-1,is)

         ! absorbing volume element into weights
         werr(:nesub,ipt,is) = werr(:nesub,ipt,is)*energy(:nesub,is)*twopi*f0_values(:nesub,is)
       end do
    end do

    ! same thing done here for lamdba as was
    ! done earlier for energy space
    do ipt=1,ng2

       if (ipt /= 1) lmodzeroes(:ipt-1) = xx(:ipt-1)
       if (ipt /= ng2) lmodzeroes(ipt:ng2-1) = xx(ipt+1:)

       call get_weights (nmax,1.0,0.0,lmodzeroes,wlerrtmp,ndiv,divmax,eflag)

       if (ipt /= 1) wlerr(:ipt-1,ipt) = wlerrtmp(:ipt-1)
       if (ipt /= ng2) wlerr(ipt+1:,ipt) = wlerrtmp(ipt:)
    end do
    deallocate(modzeroes,werrtmp,lmodzeroes,wlerrtmp)
    eflag = .false.

    !lint_error
    allocate (wlmod(-ntgrid:ntgrid,nlambda,ng2))

    do ipt = 1, ng2
       do il = 1, ng2
          wlmod(:,il,ipt) = wlerr(il,ipt)*2.0*sqrt((bmag(:)/bmax) &
               *((1.0/bmax-al(il))/(1.0/bmag(:)-al(il))))
       end do
       !If we have trapped particles use the precalculated weights
       !in wlmod as above is only for passing particles
       if (grid_has_trapped_particles()) wlmod(:,ng2+1:,ipt) = wl(:,ng2+1:)
    end do

    !How many trapped pitch angles are there?
    if (grid_has_trapped_particles() .and. new_trap_int) then
       ntrap = nlambda - ng2
       allocate(wlterr(-ntgrid:ntgrid,nlambda,ntrap))
       allocate(dummy(-ntgrid:ntgrid,nlambda))
       wlterr = 0
       call setup_trapped_lambda_grids_new_trap_int(al, dummy, wlterr)

       ! trap_error
       allocate (wtmod(-ntgrid:ntgrid,nlambda,ntrap))

       do ipt=1,ntrap
          wtmod(:,:ng2,ipt) = wl(:,:ng2)
       end do

       wtmod(:,ng2+1:,:) = wlterr(:,ng2+1:,:)
    end if

    ! eint_error
    allocate (wmod(negrid, nesub, nspec))

    wmod(:negrid-1,:,:) = werr(:,:,:)
    do is = 1,nspec
       wmod(negrid,:,is) = w(negrid,is)
    end do
  end subroutine init_weights

  !> Deallocate the weights used to provide error estimates
  subroutine finish_weights
    if (allocated(wlmod)) deallocate (wlmod)
    if (allocated(wmod)) deallocate (wmod)
    if (allocated(wtmod)) deallocate (wtmod)
  end subroutine finish_weights

  subroutine int_error_helper(g, weights, w_e, w_l, total)
    use theta_grid, only: ntgrid
    use gs2_layouts, only: g_lo, is_idx, ik_idx, it_idx, ie_idx, il_idx
    use array_utils, only: zero_array
    complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in) :: g
    real, dimension(:), intent (in) :: weights
    real, dimension(:, :), intent(in) :: w_e
    real, dimension(-ntgrid:, :), intent(in) :: w_l
    complex, dimension (-ntgrid:,:,:), intent (out) :: total
    integer :: is, il, ie, ik, it, iglo
    call zero_array(total)
    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       ik = ik_idx(g_lo,iglo)
       it = it_idx(g_lo,iglo)
       ie = ie_idx(g_lo,iglo)
       is = is_idx(g_lo,iglo)
       il = il_idx(g_lo,iglo)

       total(:, it, ik) = total(:, it, ik) + weights(is) * w_e(ie, is) * w_l(:,il) * &
            (g(:,1,iglo)+g(:,2,iglo))
    end do
  end subroutine int_error_helper

  !> FIXME : Add documentation
  subroutine lint_error (g, weights, total)
    use le_grids, only: ng2, w
    use mp, only: sum_allreduce
    implicit none
    complex, dimension(:,:,:), intent (in) :: g
    real, dimension(:), intent (in) :: weights
    complex, dimension(:,:,:,:), intent (out) :: total
    integer :: ipt
    !For each (passing) lambda point do velocity space integral
    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ipt) &
    !$OMP SHARED(ng2, g, weights, w, wlmod, total) &
    !$OMP SCHEDULE(static)
    do ipt = 1, ng2
       call int_error_helper(g, weights, w, wlmod(:, :, ipt), total(:, :, :, ipt))
    end do
    !$OMP END PARALLEL DO
    call sum_allreduce (total)
  end subroutine lint_error

  !> FIXME : Add documentation
  subroutine trap_error (g, weights, total)
    use le_grids, only: nlambda, ng2, w
    use mp, only: sum_allreduce
    implicit none
    complex, dimension(:,:,:), intent (in) :: g
    real, dimension(:), intent (in) :: weights
    complex, dimension(:,:,:,:), intent (out) :: total
    integer :: ipt
    !Loop over number of trapped points
    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ipt) &
    !$OMP SHARED(nlambda, ng2, g, weights, w, wtmod, total) &
    !$OMP SCHEDULE(static)
    do ipt = 1, nlambda - ng2
       call int_error_helper(g, weights, w, wtmod(:, :, ipt), total(:, :, :, ipt))
    end do
    !$OMP END PARALLEL DO
    call sum_allreduce (total)
  end subroutine trap_error

  !> FIXME : Add documentation
  subroutine eint_error (g, weights, total)
    use le_grids, only: nesub, wl
    use mp, only: sum_allreduce
    implicit none
    complex, dimension(:,:,:), intent (in) :: g
    real, dimension(:), intent (in) :: weights
    complex, dimension(:,:,:,:), intent (out) :: total
    integer :: ipt
    !Do velocity space integral for each ipt (for all energy grid points)
    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ipt) &
    !$OMP SHARED(nesub, g, weights, wmod, wl, total) &
    !$OMP SCHEDULE(static)
    do ipt = 1, nesub
       call int_error_helper(g, weights, wmod(:, ipt, :), wl, total(:, :, :, ipt))
    end do
    !$OMP END PARALLEL DO
    call sum_allreduce (total)
  end subroutine eint_error

  !> Error estimate based on monitoring amplitudes of legendre polynomial coefficients.
  !>
  !> FIXME: finish documentation
  subroutine get_gtran (geavg, glavg, gtavg, phi, bpar)
    use le_grids, only: nlambda, ng2, nesub, jend, grid_has_trapped_particles, il_is_trapped
    use theta_grid, only: ntgrid
    use kt_grids, only: ntheta0, naky
    use species, only: nspec
    use dist_fn_arrays, only: gnew, aj0, g_adjust, g_work, to_g_gs2, from_g_gs2
    use gs2_layouts, only: g_lo
    use array_utils, only: zero_array
    implicit none
    complex, dimension (-ntgrid:,:,:), intent (in) :: phi, bpar
    real, intent (out) :: geavg, glavg, gtavg
    real, dimension (:), allocatable :: gne2, gnl2, gnt2
    complex, dimension (:,:,:,:,:), allocatable :: getran, gltran, gttran
    real :: genorm, gemax, genum, gedenom
    real :: glnorm, glmax, glnum, gldenom
    real :: gtnorm, gtmax, gtnum, gtdenom
    integer :: ig, it, ik, is, ie, il, iglo, isgn
    allocate(gne2(0:nesub-1))
    allocate(gnl2(0:ng2-1))

    genorm = 0.0 ; glnorm = 0.0
    gne2  = 0.0 ; gnl2 = 0.0
    gemax = 0.0 ; glmax = 0.0
    genum = 0.0 ; gedenom = 0.0
    glnum = 0.0 ; gldenom = 0.0
    gtnum = 0.0 ; gtdenom = 0.0

    allocate(getran(0:nesub-1, -ntgrid:ntgrid, ntheta0, naky, nspec))
    allocate(gltran(0:ng2-1, -ntgrid:ntgrid, ntheta0, naky, nspec))

    call zero_array(getran); call zero_array(gltran)

    ! Only needed if we have trapped particles, but always allocate to support
    ! omp reduction later etc.
    allocate(gnt2(0:2*(nlambda-ng2-1)))
    allocate(gttran(0:2*(nlambda-ng2-1), -ntgrid:ntgrid, ntheta0, naky, nspec))
    gtnorm = 0.0 ; gnt2 = 0.0 ; gtmax = 0.0 ; call zero_array(gttran)

    ! transform from g to h
    call g_adjust (gnew, phi, bpar, direction = from_g_gs2)

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(iglo, isgn, ig) &
    !$OMP SHARED(g_lo, ntgrid, g_work, aj0, gnew) &
    !$OMP COLLAPSE(3) &
    !$OMP SCHEDULE(static)
    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       do isgn = 1, 2
          do ig = -ntgrid, ntgrid
             g_work(ig, isgn, iglo) = aj0(ig, iglo) * gnew(ig, isgn, iglo)
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    ! transform from h back to g
    call g_adjust (gnew, phi, bpar, direction = to_g_gs2)

    ! perform legendre transform on dist. fn. to obtain coefficients
    ! used when expanding dist. fn. in legendre polynomials
    if (grid_has_trapped_particles()) then
       call legendre_transform (g_work, getran, gltran, gttran)
    else
       call legendre_transform (g_work, getran, gltran)
    end if

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(is, ik, it, ig, ie, il, gne2, gnl2, genorm, gemax, glnorm, glmax, &
    !$OMP gnt2, gtnorm, gtmax) &
    !$OMP SHARED(nspec, naky, ntheta0, ntgrid, getran, gltran, gttran, nesub, ng2, jend) &
    !$OMP REDUCTION(+: genum, gedenom, glnum, gldenom, gtnum, gtdenom) &
    !$OMP SCHEDULE(static)
    do is = 1, nspec
       do ik = 1, naky
          do it = 1, ntheta0
             do ig = -ntgrid, ntgrid
                do ie = 0, nesub-1
                   gne2(ie) = real(getran(ie,ig,it,ik,is)*conjg(getran(ie,ig,it,ik,is)))
                end do

                do il=0,ng2-1
                   gnl2(il) = real(gltran(il,ig,it,ik,is)*conjg(gltran(il,ig,it,ik,is)))
                end do

                genorm = maxval(gne2)
                if (nesub < 3) then
                   gemax = gne2(size(gne2)-1)
                else
                   gemax = maxval(gne2(nesub-3:nesub-1))
                end if
                glnorm = maxval(gnl2)
                glmax = maxval(gnl2(ng2-3:ng2-1))

                genum = genum + gemax
                gedenom = gedenom + genorm
                glnum = glnum + glmax
                gldenom = gldenom + glnorm

                if (grid_has_trapped_particles()) then
                   do il=0, 2*(jend(ig)-ng2-1)
                      gnt2(il) = real(gttran(il,ig,it,ik,is)*conjg(gttran(il,ig,it,ik,is)))
                   end do
                   gtnorm = maxval(gnt2(0:2*(jend(ig)-ng2-1)))
                   if (il_is_trapped(jend(ig))) then
                      gtmax = maxval(gnt2(2*(jend(ig)-ng2-1)-2:2*(jend(ig)-ng2-1)))
                   else
                      gtmax = gnt2(0)
                   end if

                   gtnum = gtnum + gtmax
                   gtdenom = gtdenom + gtnorm
                end if
             end do
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    geavg = genum / gedenom
    glavg = glnum / gldenom
    if (grid_has_trapped_particles()) gtavg = gtnum / gtdenom

    deallocate(gne2, gnl2, getran, gltran, gnt2, gttran)
  end subroutine get_gtran

  subroutine setup_legendre_transform
    use le_grids, only: nesub, nlambda, ng2, vcut, speed, jend, al, negrid, xx
    use le_grids, only: grid_has_trapped_particles
    use theta_grid, only: bmag, bmax, ntgrid
    use species, only: nspec
    real, dimension (:), allocatable :: nodes
    real, dimension (:,:), allocatable :: lpltmp, lpttmp
    real :: ulim
    integer :: is, il, ig, ntrap, lpesize
    lpesize = nesub

    if (allocated(lpl)) return
    allocate(lpltmp(ng2,0:ng2-1))
    allocate(lpl(nlambda,0:ng2-1))
    allocate(lpe(negrid,0:lpesize-1,nspec)) ; lpe = 0.0

    ! get value of first nesub legendre polynomials
    ! at each of the grid points on (0,vcut)
    do is = 1, nspec
       call legendre_polynomials (0.0,vcut,speed(:lpesize,is),lpe(:lpesize,:,is))
    end do

    ! get value of first ng2 legendre polynomials
    ! at each of the grid points on (0,1)
    call legendre_polynomials (0.0,1.0,xx,lpltmp)

    lpl = 0.0
    lpl(1:ng2,:) = lpltmp

    if (grid_has_trapped_particles()) then
       allocate (lpt(nlambda,0:2*(nlambda-ng2-1),-ntgrid:ntgrid,2))
       lpt = 0.0
       do ig = -ntgrid, ntgrid
          ntrap = 1
          if (jend(ig) > ng2+1) then
             ntrap = jend(ig)-ng2
             allocate (nodes(2*ntrap-1))
             allocate (lpttmp(2*ntrap-1,0:2*(ntrap-1)))
             do il = 1, ntrap
                nodes(il) = -sqrt(max(0.0,1.0-al(ng2+il)*bmag(ig)))
             end do
             nodes(ntrap+1:) = -nodes(ntrap-1:1:-1)

             ulim = sqrt(1.0-bmag(ig)/bmax)
             call legendre_polynomials (-ulim,ulim,nodes,lpttmp)
             lpt(ng2+1:jend(ig),0:2*(ntrap-1),ig,2) = lpttmp(1:ntrap,:)
             lpt(ng2+1:jend(ig)-1,0:2*(ntrap-1),ig,1) = lpttmp(2*ntrap-1:ntrap+1:-1,:)
             deallocate (nodes, lpttmp)
          end if
       end do
    end if

    deallocate (lpltmp)

  end subroutine setup_legendre_transform

  subroutine finish_legendre_transform
    if (allocated(lpl)) deallocate(lpl, lpe, lpt)
  end subroutine finish_legendre_transform
  
  !> FIXME : Add documentation
  !>
  !> @note This routine depends on the values stored in [[xx]] as input
  !> to the legendre_polynomials method. However, with Radau grids the
  !> values in [[xx]] are not Legendre roots so this routine is probably
  !> incorrect when [[le_grids_knobs::radau_gauss_grid]] is `.true.`.
  subroutine legendre_transform (g, tote, totl, tott)  
    use mp, only: sum_allreduce
    use theta_grid, only: ntgrid
    use gs2_layouts, only: g_lo, it_idx, ik_idx, il_idx, ie_idx, is_idx
    use le_grids, only: nesub, w, wl, ng2, jend
    use array_utils, only: zero_array
    implicit none
    complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in) :: g
    complex, dimension (0:,-ntgrid:,:,:,:), intent (out) :: tote, totl
    complex, dimension (0:,-ntgrid:,:,:,:), intent (out), optional :: tott
    complex :: totfac
    integer :: is, il, ie, ik, it, iglo, ig, im, lpesize

    lpesize = nesub
    ! carry out legendre transform to get coefficients of
    ! legendre polynomial expansion of g
    call zero_array(tote) ; call zero_array(totl)

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(iglo, it, ik, il, ie, is, totfac, im) &
    !$OMP SHARED(g_lo, ntgrid, w, wl, g, lpesize, ng2, lpe, lpl) &
    !$OMP REDUCTION(+: tote, totl) &
    !$OMP SCHEDULE(static)
    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       it = it_idx(g_lo, iglo) ; ik = ik_idx(g_lo, iglo)
       il = il_idx(g_lo, iglo) ; ie = ie_idx(g_lo, iglo) ; is = is_idx(g_lo, iglo)
       do ig = -ntgrid, ntgrid
          totfac = w(ie,is) * wl(ig,il) * (g(ig,1,iglo) + g(ig,2,iglo))
          do im = 0, lpesize - 1
             tote(im, ig, it, ik, is) = tote(im, ig, it, ik, is) &
                  + totfac * lpe(ie,im,is) * (2*im+1)
          end do
          do im = 0, ng2 - 1
             totl(im, ig, it, ik, is) = totl(im, ig, it, ik, is) &
                  + totfac * lpl(il,im) * (2*im+1)
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !Now complete velocity integral, bringing back results to proc0
    call sum_allreduce (tote)
    call sum_allreduce (totl)

    if (present(tott)) then
       call zero_array(tott)
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(iglo, it, ik, il, ie, is, im) &
       !$OMP SHARED(g_lo, ntgrid, w, wl, g, jend, ng2, lpt) &
       !$OMP REDUCTION(+: tott) &
       !$OMP SCHEDULE(static)
       do iglo = g_lo%llim_proc, g_lo%ulim_proc
          it = it_idx(g_lo, iglo) ; ik = ik_idx(g_lo, iglo)
          il = il_idx(g_lo, iglo) ; ie = ie_idx(g_lo, iglo) ; is = is_idx(g_lo, iglo)
          do ig = -ntgrid, ntgrid
             do im = 0, 2 * (jend(ig)-ng2-1)
                tott(im, ig, it, ik, is) = tott(im, ig, it, ik, is) + &
                     w(ie,is) * wl(ig,il) * (lpt(il,im,ig,1) * g(ig,1,iglo) &
                     + lpt(il,im,ig,2) * g(ig,2,iglo)) * (2*im+1)
             end do
          end do
       end do
       !$OMP END PARALLEL DO
       call sum_allreduce (tott)
    end if
  end subroutine legendre_transform

  !> FIXME : Add documentation  
  subroutine legendre_polynomials (llim, ulim, xptsdum, lpdum)
    use iso_fortran_env, only: real64
    implicit none
    real, intent (in) :: ulim, llim
    real, dimension (:), intent (in)   :: xptsdum
    real, dimension (:,0:), intent(out) :: lpdum
    real(real64), dimension (:), allocatable :: lp1, lp2, lp3, zshift
    integer :: j, mmax

    lpdum = 0.0
    mmax = size(xptsdum)

    allocate(lp1(mmax),lp2(mmax),lp3(mmax),zshift(mmax))

    lp1 = real(1.0,kind(lp1(1)))
    lp2 = real(0.0,kind(lp2(1)))

    lpdum(:,0) = real(1.0,kind(lpdum))

    zshift = real(2.0,kind(zshift))*(xptsdum-llim)/(ulim-llim) - real(1.0,kind(zshift))

    do j=1, size(lpdum(1,:))-1
       lp3 = lp2
       lp2 = lp1
       lp1 = ((2*j-1) * zshift * lp2 - (j-1) * lp3) / j
       lpdum(:,j) = lp1
    end do

    deallocate(lp1,lp2,lp3,zshift)
  end subroutine legendre_polynomials
end module diagnostics_velocity_space

gyrokinetics / gs2 / 1970319704

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous