init_vpdiff Subroutine

  subroutine init_vpdiff
    use le_grids, only: al, nlambda, jend, ng2, il_is_wfb, il_is_passing, ixi_to_il, ixi_to_isgn
    use theta_grid, only: ntgrid, bmag
    use array_utils, only: zero_array
    implicit none

    integer :: il, isgn, ixi, ig, je, te
    real :: slb0, slb1, slb2, slbl, slbr

    if (.not. allocated(vpdiff) .and. conservative) then

       allocate (vpdiff(-ntgrid:ntgrid, 2, nlambda)) ; call zero_array(vpdiff)
       
       do ig = -ntgrid, ntgrid
          je = jend(ig)
          if (il_is_passing(je) .or. (il_is_wfb(je) .and. special_wfb_lorentz)) then
             te = ng2
          else
             te = je
          end if
          do il = 2, te-1
             slb0 = safe_sqrt(1.0 - bmag(ig) * al(il-1))
             slb1 = safe_sqrt(1.0 - bmag(ig) * al(il))
             slb2 = safe_sqrt(1.0 - bmag(ig) * al(il+1))
             
             slbl = (slb1 + slb0) * 0.5  ! xi(j-1/2)
             slbr = (slb1 + slb2) * 0.5  ! xi(j+1/2)
             
             vpdiff(ig, 1, il) = (slbl**2 - slbr**2) / pitch_weights(il, ig)
          end do

          ! boundary at xi = 1
          slb1 = safe_sqrt(1.0 - bmag(ig) * al(1))
          slb2 = safe_sqrt(1.0 - bmag(ig) * al(2))
          slbr = 0.5 * (slb1 + slb2)
          vpdiff(ig, 1, 1) = (1.0 - slbr**2) / pitch_weights(1, ig)
          
          ! boundary at xi = 0
          il = te
          slb0 = safe_sqrt(1.0 - bmag(ig) * al(il-1))
          if (te == ng2) then ! Passing
             slb1 = safe_sqrt(1.0 - bmag(ig) * al(il))
             slb2 = -slb1
          else
             ! We would expect safe_sqrt(1.0 - bmag(ig) * al(il)) = 0 here so could just
             ! use slb1 = safe_sqrt(1.0 - bmag(ig) * al(il)) = 0 for both branches?
             slb1 = 0.0
             slb2 = -slb0
          end if
          
          slbl = (slb1 + slb0) * 0.5
          slbr = (slb1 + slb2) * 0.5

          if (abs(pitch_weights(il, ig)) <= epsilon(0.0)) then
             vpdiff(ig, 1, il) = 0.0
          else
             vpdiff(ig, 1, il) = (slbl**2 - slbr**2) / pitch_weights(il, ig)
          end if
          
          vpdiff(ig, 2, :) = -vpdiff(ig, 1, :)
          
       end do

       if (use_le_layout) then
          allocate(vpdiffle(nxi_lim, -ntgrid:ntgrid))

          do ig = -ntgrid, ntgrid
             do ixi = 1, nxi_lim
                il = ixi_to_il(ixi, ig)
                isgn = ixi_to_isgn(ixi, ig)
                vpdiffle(ixi, ig) = vpdiff(ig, isgn, il)
             end do
          end do
       end if
    end if
    
  end subroutine init_vpdiff