hahm_burrell Subroutine

  subroutine hahm_burrell(self, a)
    use file_utils, only: open_output_file, close_output_file
    use mp, only: proc0
    implicit none
    class(abstract_geo_type), intent(in) :: self
    real, intent(in) :: a
    integer :: i, fort24_unit, nr
    real :: gradpsi, mag_B, rho_eq, rp1, rp2, rho1, rho2, drhodpsiq
    real, dimension(self%nr) :: gamma_shear, pbar, dp, d2p, pres, bs_coll, s__hat, bs_sh
    if (.not. proc0) return
    nr = self%nr
    gamma_shear = 0.
    pbar = (self%eqpsi - self%eqpsi(1)) / (self%eqpsi(nr) - self%eqpsi(1))

    do i = 1, nr
       dp(i)  = self%dpfun(pbar(i))
       pres(i) = self%pfun(pbar(i))
    end do

    d2p(1) = 0.0 ; d2p(nr) = 0.0
    do i = 2, nr - 1
       d2p(i) = (dp(i+1) - dp(i-1)) / (self%eqpsi(i+1) - self%eqpsi(i-1))
    end do

    do i = 2, nr- 1
       rp1 = self%eqpsi(i+1) ; rp2 = self%eqpsi(i-1)
       rho1 = 0.5 * self%diameter(rp1) ; rho2 = 0.5 * self%diameter(rp2)
       drhodpsiq = (rho1 - rho2) / (rp1 - rp2)
       gradpsi = self%eqitem(self%eqpsi(i), 0., self%dpbish(:,:,1), 'R')

       mag_B = hypot(self%btori(pbar(i)), gradpsi) * self%invR(self%eqpsi(i), 0.)
       !For geq and ideq could be mag_B = geom%eqitem(geom%eqpsi(i), 0., geom%B_psi, 'R')

       gamma_shear(i) = 0.01 * gradpsi**2 * (d2p(i) / pres(i) - a * (dp(i) / pres(i))**2) &
            * (2 * pres(i) / self%beta_0)**((1 - a) / 2) * (-pres(i) / (dp(i) / drhodpsiq)) &
            / mag_B

       s__hat(i) = (self%qfun(pbar(i+1)) - self%qfun(pbar(i-1))) / (rho2 - rho1)
       s__hat(i) = s__hat(i) * 0.5 * (rho1 + rho2) / self%qfun(pbar(i))
    end do

    ! Assume rho_*0 = 0.01, nu_*0 = 1.0
    do i = 2, nr - 1
       bs_coll(i) = 0.01 * ((2 * pres(i) / self%beta_0)**(1 - a))**2.5 &
             * self%qfun(0.) * 0.5 * self%diameter(self%eqpsi(i)) &
            * (-dp(i) / pres(i)) * self%R_psi(1,1)**1.5 / (2 * pres(i) / self%beta_0)**a
       bs_sh(i) = (-pres(i) / (dp(i) / drhodpsiq))*s__hat(i) / (self%R_psi(1,1) &
            * self%qfun(pbar(i)) * sqrt(pres(i)))
    end do

    call open_output_file(fort24_unit,".fort.24")
    do i = 2, nr - 1
       rho_eq = 0.5 * self%diameter(self%eqpsi(i))
       write(fort24_unit,'(i5,11(1x,e16.9))') i, pbar(i), rho_eq, pres(i), gamma_shear(i), bs_coll(i), bs_sh(i)
    end do
    call close_output_file(fort24_unit)
  end subroutine hahm_burrell