drift Subroutine

  subroutine drift(rgrid, rp, gradstot, gradrptot, gradztot, gradbtot, dqdrp, dpsidrho, &
       dpsidrp, bmag, invR, theta, dbetadrho, bi, force_sym, nperiod, &
       gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0)
    implicit none
    real, dimension (-ntgrid:), intent (in) :: rgrid, bmag, invR, theta
    real, dimension (-ntgrid:, :), intent (in) :: gradstot, gradrptot, gradztot, gradbtot
    real, intent (in) :: rp, dqdrp, dpsidrho, dpsidrp, dbetadrho, bi
    logical, intent(in) :: force_sym
    integer, intent(in) :: nperiod
    real, dimension(-ntgrid:ntgrid), intent(out) :: gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0
    real, dimension (-ntgrid:ntgrid, 3) :: pgradtot, dummy, curve, bvector
    real, dimension (-ntgrid:ntgrid) :: gbdrift1, gbdrift2, cvdrift1, cvdrift2, cdrift1, cdrift2
    real :: dp
    if (use_large_aspect) then
       gbdrift = (-sin(theta) * gradstot(:,1) - cos(theta) * gradstot(:,2)) * invR
       gbdrift = 2 * dpsidrho * gbdrift

       gbdrift0 = (-sin(theta) * gradrptot(:,1) - cos(theta) * gradrptot(:,2)) * invR
       gbdrift0 = 2 * dpsidrho * gbdrift0 * dqdrp

       cvdrift = gbdrift ; cvdrift0 = gbdrift0
    else
       !     compute magnetic field
       call bvectortgrid(invR, nth, gradrptot, dpsidrp, bi, bvector)
       dp = dbetadrho / (2 * dpsidrho)

       call crosstgrid(bvector, gradbtot, dummy) !B cross Grad B
       call periodic_copy(dummy(:, 1), 0.0, nperiod)
       call periodic_copy(dummy(:, 2), 0.0, nperiod)
       call periodic_copy(dummy(:, 3), 0.0, nperiod)
       call dottgridf(dummy, gradstot, gbdrift1)
       call dottgridf(dummy, gradrptot, gbdrift2)

       !     create curvature
       if (bishop == 0) then
          call geom%gradient(rgrid, theta, pgradtot(:, 1:2), 'R', rp, nth, ntgrid, bishop /= 0)
          pgradtot(-nth:nth, 3) = 0.0
       else
          pgradtot = gradrptot * dpsidrp * dp !Grad rho * dp/drho
       end if

       curve = gradbtot
       curve(-nth:nth, 1) = curve(-nth:nth, 1) + pgradtot(-nth:nth, 1) / bmag(-nth:nth)
       curve(-nth:nth, 2) = curve(-nth:nth, 2) + pgradtot(-nth:nth, 2) / bmag(-nth:nth)

       call crosstgrid(bvector, curve, dummy) !B cross [Grad B + Grad P / B]
       call periodic_copy(dummy(:, 1), 0.0, nperiod)
       call periodic_copy(dummy(:, 2), 0.0, nperiod)
       call periodic_copy(dummy(:, 3), 0.0, nperiod)
       call dottgridf(dummy, gradstot, cvdrift1)
       call dottgridf(dummy, gradrptot, cvdrift2)

       ! factor of 2 appears below because will later multiply
       ! by wunits, which has a factor of 1/2 built-in
       gbdrift = 2 * dpsidrho * gbdrift1 / bmag**3
       gbdrift0 = 2 * dpsidrho * gbdrift2 * dqdrp / bmag**3
       cvdrift = 2 * dpsidrho * cvdrift1 / bmag**3
       cvdrift0 = 2 * dpsidrho * cvdrift2 * dqdrp / bmag**3
    end if

    ! coriolis drift term is -2*vpa*om_phi/Omega * (bhat x (zhat x bhat)) . grad (g+q<phi>F_0/T)
    ! here we calculate 2*(zhat . grad (alpha+q*theta0))/(B/Bref)

    ! this is the Zhat . grad(alpha) part of Zhat . grad(S)
    ! note that a factor of k_y = (n0/a) (drho/dpsiN) has been factored out
    call dottgridf(gradztot, gradstot, cdrift1)
    ! this is the Zhat . grad(q) part of Zhat . grad(S)
    ! note that a factor of k_y * theta0 has been factored out
    call dottgridf(gradztot, gradrptot, cdrift2)

    ! factor of 4 (instead of 2) below because will later mutliply
    ! by wunits, which has a built-in factor of 1/2
    ! cdrift0 is part of coriolis drift dependent on theta0
    cdrift  = -4 * dpsidrho * cdrift1 / bmag
    cdrift0 = -4 * dpsidrho * cdrift2 * dqdrp / bmag

    if (force_sym) then
       call sym(gbdrift, 0, ntgrid)
       call sym(cvdrift, 0, ntgrid)
       call sym(gbdrift0, 1, ntgrid)
       call sym(cvdrift0, 1, ntgrid)
       ! should maybe add in cdrift and cdrift0 here -- MAB
    end if
  end subroutine drift