initialize_GC_fields Subroutine

  subroutine initialize_GC_fields(F)
    !! Computes the auxiliary fields \(\nabla|{\bf B}|\) and
    !! \(\nabla\times\hat{b}\) that are used in the RHS of the
    !! evolution equations for the GC orbit model.
    TYPE(FIELDS), INTENT(INOUT)      :: F
    !! An instance of the KORC derived type FIELDS.
    INTEGER                        :: ii
    !! Iterator across F%dim
    REAL(rp), DIMENSION(:,:),ALLOCATABLE :: Bmag
    !! Magnetic field magnitude
    REAL(rp), DIMENSION(:,:,:),ALLOCATABLE :: bhat
    !! Magnetic field unit vector

    Bmag=SQRT(F%B_2D%R**2+F%B_2D%PHI**2+F%B_2D%Z**2)

    ALLOCATE(bhat(F%dims(1),F%dims(3),3))

    bhat(:,:,1)=F%B_2D%R/Bmag
    bhat(:,:,2)=F%B_2D%PHI/Bmag
    bhat(:,:,3)=F%B_2D%Z/Bmag


    F%gradB_2D%PHI=0.
    ! No variation in phi direction

    ! Single-sided difference for axiliary fields at edge nodes
    ! Differential over R on first index, differential over Z
    ! on second index.

    ! gradB
    ! edge nodes at minimum R,Z
    F%gradB_2D%R(1,:)=(Bmag(2,:)-Bmag(1,:))/(F%X%R(2)-F%X%R(1))
    F%gradB_2D%Z(:,1)=(Bmag(:,2)-Bmag(:,1))/(F%X%Z(2)-F%X%Z(1))

    ! edge nodes at maximum R,Z
    F%gradB_2D%R(F%dims(1),:)=(Bmag(F%dims(1),:)-Bmag(F%dims(1)-1,:))/ &
         (F%X%R(F%dims(1))-F%X%R(F%dims(1)-1))
    F%gradB_2D%Z(:,F%dims(3))=(Bmag(:,F%dims(3))-Bmag(:,F%dims(3)-1))/ &
         (F%X%Z(F%dims(3))-F%X%Z(F%dims(3)-1))

    ! curlb
    ! edge nodes at minimum R,Z
    ! R component has differential over Z
    F%curlb_2D%R(:,1)=-(bhat(:,2,2)-bhat(:,1,2))/ &
         (F%X%Z(2)-F%X%Z(1))

    ! PHI component has differentials over R and Z
    F%curlb_2D%PHI(1,:)=-(bhat(2,:,3)-bhat(1,:,3))/ &
         (F%X%R(2)-F%X%R(1))         

    F%curlb_2D%PHI(:,1)=F%curlb_2D%PHI(:,1)+ &
         ((bhat(:,2,1)-bhat(:,1,1))/(F%X%Z(2)-F%X%Z(1)))

    ! Z component has differentials over R
    F%curlb_2D%Z(1,:)=((bhat(2,:,2)*F%X%R(2)- &
         bhat(1,:,2)*F%X%R(1))/(F%X%R(2)-F%X%R(1)))/F%X%R(1)

    ! edge nodes at minimum R,Z
    ! R component has differential over Z
    F%curlb_2D%R(:,F%dims(3))=-(bhat(:,F%dims(3),2)- &
         bhat(:,F%dims(3)-1,2))/ &
         (F%X%Z(F%dims(3))-F%X%Z(F%dims(3)-1))

    ! PHI component has differentials over R and Z
    F%curlb_2D%PHI(F%dims(1),:)=F%curlb_2D%PHI(F%dims(1),:)- &
         (bhat(F%dims(1),:,3)-bhat(F%dims(1)-1,:,3))/ &
         (F%X%R(F%dims(1))-F%X%R(F%dims(1)-1))         

    F%curlb_2D%PHI(:,F%dims(3))=F%curlb_2D%PHI(:,F%dims(3))+ &
         ((bhat(:,F%dims(3),1)-bhat(:,F%dims(3)-1,1))/ &
         (F%X%Z(F%dims(3))-F%X%Z(F%dims(3)-1)))

    ! Z component has differentials over R
    F%curlb_2D%Z(F%dims(1),:)=((bhat(F%dims(1),:,2)*F%X%R(F%dims(1))- &
         bhat(F%dims(1)-1,:,2)*F%X%R(F%dims(1)-1))/(F%X%R(F%dims(1))- &
         F%X%R(F%dims(1)-1)))/F%X%R(F%dims(1))

    do ii=2_idef,F%dims(1)-1
       ! central difference over R for interior nodes
       F%gradB_2D%R(ii,:)=(Bmag(ii+1,:)-Bmag(ii-1,:))/ &
            (F%X%R(ii+1)-F%X%R(ii-1))
       F%curlb_2D%Z(ii,:)=((bhat(ii+1,:,2)*F%X%R(ii+1)- &
            bhat(ii-1,:,2)*F%X%R(ii-1))/(F%X%R(ii+1)-F%X%R(ii-1)))/ &
            F%X%R(ii)
       F%curlb_2D%PHI(ii,:)=F%curlb_2D%PHI(ii,:)- &
            (bhat(ii+1,:,3)-bhat(ii-1,:,3))/ &
            (F%X%R(ii+1)-F%X%R(ii-1))
    end do
    do ii=2_idef,F%dims(3)-1
       ! central difference over Z for interior nodes
       F%gradB_2D%Z(:,ii)=(Bmag(:,ii+1)-Bmag(:,ii-1))/ &
            (F%X%Z(ii+1)-F%X%Z(ii-1))
       F%curlb_2D%R(:,ii)=-(bhat(:,ii+1,2)-bhat(:,ii-1,2))/ &
            (F%X%Z(ii+1)-F%X%Z(ii-1))
       F%curlb_2D%PHI(:,ii)=F%curlb_2D%PHI(:,ii)+ &
            ((bhat(:,ii+1,1)-bhat(:,ii-1,1))/(F%X%Z(ii+1)-F%X%Z(ii-1)))
    end do

    DEALLOCATE(Bmag)
    DEALLOCATE(bhat) 

  end subroutine initialize_GC_fields