pure subroutine calculate_metric_terms(gradstot, gradrptot, gradthtot, dpsidrho, dqdrp, &
qval, rhoc, bmag, force_sym, gds2, gds21, gds22, gds23, gds24, gds24_noq)
real, dimension(-ntgrid:, :), intent(in) :: gradstot, gradrptot, gradthtot
real, intent(in) :: dpsidrho, dqdrp, qval, rhoc
real, dimension(-ntgrid:), intent(in) :: bmag
logical, intent(in) :: force_sym
real, dimension(-ntgrid:), intent(out) :: gds2, gds21, gds22, gds23, gds24, gds24_noq
real, dimension(:), allocatable :: gdsdum1, gdsdum2, gdsdum3, gdsdum4, gdsdum5
allocate(gdsdum1(-ntgrid:ntgrid), gdsdum2(-ntgrid:ntgrid), gdsdum3(-ntgrid:ntgrid))
call dottgridf(gradstot, gradstot, gdsdum1)
call dottgridf(gradstot, gradrptot, gdsdum2)
call dottgridf(gradrptot, gradrptot, gdsdum3)
gds2 = gdsdum1 *dpsidrho**2
gds21 = gdsdum2*dqdrp *dpsidrho**2
gds22 = gdsdum3*dqdrp*dqdrp*dpsidrho**2
! v_E . grad theta = i*ky*phi*(rho_{r}*v_{t,r}/a**2)*(gds23 + gsd24*theta0)
! calculation of gds23 and gds24 assumes rp=rho, which is not necessarily true
! for numerical equilibria
allocate(gdsdum4(-ntgrid:ntgrid), gdsdum5(-ntgrid:ntgrid))
call dottgridf(gradrptot, gradthtot, gdsdum4) !grad rho dot grad theta
call dottgridf(gradstot, gradthtot, gdsdum5) !grad alpha dot grad theta
gds23 = (gdsdum4 * gds2 - gdsdum5 * gdsdum2 * dpsidrho**2) / bmag**2
gds24_noq = (gdsdum4 * gdsdum2 - gdsdum5 * gdsdum3) * (dpsidrho / bmag)**2
gds24 = gds24_noq * dqdrp
gds24_noq = gds24_noq * qval / rhoc
if (force_sym) then
call sym(gds2, 0, ntgrid)
call sym(gds21, 1, ntgrid)
call sym(gds22, 0, ntgrid)
! not sure about the following yet -- MAB
! call sym(gds23, 0, ntgrid); call sym(gds24, 1, ntgrid)
end if
end subroutine calculate_metric_terms