calculate_surface_area Subroutine

  subroutine calculate_surface_area(theta, gradpar, bmag, grho, dpsidrho, rhoc, nth, &
       surfarea, dvdrhon)
    use constants, only: twopi
    use file_utils, only: open_output_file, close_output_file
    real, dimension(-ntgrid:), intent(in) :: theta, gradpar, bmag, grho
    real, intent(in) :: dpsidrho, rhoc
    integer, intent(in) :: nth
    real, intent(out) :: surfarea, dvdrhon
    real, dimension(:), allocatable :: jacob, ds, dummy
    integer :: i, fort25_unit
    ! compute grad rho*Jacobian of the transformation to the new coordinates:
    ! Axisymmetry assumed to get this (nu indep. of phi):
    allocate(jacob(-ntgrid:ntgrid), ds(-ntgrid:ntgrid), dummy(-ntgrid:ntgrid))
    do i = -nth, nth
       jacob(i) = abs(dpsidrho / (gradpar(i)*bmag(i)))
       ds(i) = grho(i) * jacob(i)
    end do

    ! compute surface area from A=Int(R sqrt(r**2+drdtheta**2) dtheta dphi)
    !  does not work with equal_arc grid
    ! compute the surface area from A=Int(J |grad rho| dtheta dalpha)
    call integrate(ds, theta, dummy, nth)
    surfarea = twopi * (dummy(nth) - dummy(-nth))

    ! compute dV/drhon = 2 pi Int(J dtheta)
    call integrate(jacob, theta, dummy, nth)
    dvdrhon = twopi * (dummy(nth) - dummy(-nth))

    call open_output_file(fort25_unit,".fort.25")
    write(fort25_unit,*) rhoc, f_trap(bmag(-nth:nth), theta(-nth:nth), jacob(-nth:nth))
    call close_output_file(fort25_unit)
  end subroutine calculate_surface_area