setup_get_fields_direct_from_dfn_denominators Subroutine

  subroutine setup_get_fields_direct_from_dfn_denominators(phi_denom, bpar_denom)
    use theta_grid, only: ntgrid, bmag
    use kt_grids, only: naky, ntheta0
    use run_parameters, only: beta, has_phi, has_bpar
    use array_utils, only: copy
    implicit none
    real, dimension(:, :, :), allocatable, intent(in out) :: phi_denom, bpar_denom
    integer :: it, ik, ig

    ! get phi denominator
    if (has_phi) then
       ! Note we allocate and set as the full size to support calls to
       ! this routine with either choice for local_only.
       if (.not. allocated(phi_denom)) allocate(phi_denom(-ntgrid:ntgrid, ntheta0, naky))

       if (has_bpar) then
          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(it, ik) &
          !$OMP SHARED(ntheta0, naky, beta, gamtot2, &
          !$OMP bmag, gamtot1, phi_denom, gamtot) &
          !$OMP COLLAPSE(2) &
          !$OMP SCHEDULE(static)
          do ik = 1, naky
             do it = 1, ntheta0
                phi_denom(:, it, ik) = (beta * gamtot2(:, it, ik) + bmag**2) &
                     * gamtot(:, it, ik) &
                     + (beta/2.0) * gamtot1(:, it, ik) * gamtot1(:, it, ik)
             end do
          end do
          !$OMP END PARALLEL DO
       else
          ! If fbpar = 0 then the terms in earlier branch proportional to beta
          ! should be zero as these terms arrive from the bpar parts of quasineutrality
          ! as expressed in terms of gnew. In this case we get the below.
          call copy(gamtot, phi_denom)
       end if

       ! Now invert the denominator
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(it, ik, ig) &
       !$OMP SHARED(ntheta0, naky, ntgrid, phi_denom) &
       !$OMP COLLAPSE(2) &
       !$OMP SCHEDULE(static)
       do ik = 1, naky
          do it = 1, ntheta0
             do ig = -ntgrid, ntgrid
                if (abs(phi_denom(ig, it, ik)) > epsilon(0.0)) then
                   phi_denom(ig, it, ik) = &
                        1.0 / phi_denom(ig, it, ik)
                else
                   phi_denom(ig, it, ik) = 0.0
                end if
             end do
          end do
       end do
       !$OMP END PARALLEL DO

    end if

    ! get bpar
    if (has_bpar) then
       ! Note we allocate and set as the full size to support calls to
       ! this routine with either choice for local_only.
       if (.not. allocated(bpar_denom)) allocate(bpar_denom(-ntgrid:ntgrid, ntheta0, naky))

       ! As in our calculation of phi, we see that the equations for
       ! bpar and phi are coupled when written in terms of gnew. This
       ! means we get contributions arising from phi in the expression
       ! for B|| here. In particular, they are the terms involving
       ! gamtot1. However, if we have disabled phi then they should
       ! not appear in our calculation, so handle that here.
       if (has_phi) then
          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(it, ik) &
          !$OMP SHARED(ntheta0, naky, beta, gamtot, &
          !$OMP gamtot1, bpar_denom, gamtot2, bmag) &
          !$OMP COLLAPSE(2) &
          !$OMP SCHEDULE(static)
          do ik = 1, naky
             do it = 1, ntheta0
                bpar_denom(:, it, ik) = gamtot(:, it, ik) * &
                     (beta * gamtot2(:, it, ik) + bmag**2) + &
                     (beta/2.0) * gamtot1(:, it, ik) * gamtot1(:, it, ik)
             end do
          end do
          !$OMP END PARALLEL DO
       else
          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(it, ik) &
          !$OMP SHARED(ntheta0, naky, beta, &
          !$OMP bpar_denom, gamtot2, bmag) &
          !$OMP COLLAPSE(2) &
          !$OMP SCHEDULE(static)
          do ik = 1, naky
             do it = 1, ntheta0
                bpar_denom(:, it, ik) = (beta * gamtot2(:, it, ik) + bmag**2)
             end do
          end do
          !$OMP END PARALLEL DO
       end if

       ! Now invert the denominator
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(it, ik, ig) &
       !$OMP SHARED(ntheta0, naky, ntgrid, bpar_denom) &
       !$OMP COLLAPSE(2) &
       !$OMP SCHEDULE(static)
       do ik = 1, naky
          do it = 1, ntheta0
             do ig = -ntgrid, ntgrid
                if (abs(bpar_denom(ig, it, ik)) > epsilon(0.0)) then
                   bpar_denom(ig, it, ik) = &
                        1.0 / bpar_denom(ig, it, ik)
                else
                   bpar_denom(ig, it, ik) = 0.0
                end if
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    end if
  end subroutine setup_get_fields_direct_from_dfn_denominators