check_theta_grid_salpha Subroutine

  subroutine check_theta_grid_salpha(report_unit,alne,dbetadrho)
    use theta_grid_params, only: eps, epsl, pk, shat
    use warning_helpers, only: exactly_equal
    implicit none
    integer, intent(in) :: report_unit
    real, intent(in) :: alne, dbetadrho
!CMR input dbetadrho is computed externally (eg from species)
!    allowing consistency check
    real :: arat, qsf
!
 ! Find q, r/R, R/a
 !
    if (epsl > 0.) then
       arat = 2. / epsl

       if (exactly_equal(epsl, 2.0)) then
          write (report_unit, &
               & fmt="('Scale lengths are normalized to the major radius, R')")
       else
          write (report_unit, fmt="('The aspect ratio R/a = ',f7.4)") arat
          if (exactly_equal(alne, 1.0)) then
             write (report_unit, &
                  & fmt="('Scale lengths are normalized to the density scale length, Ln')")
          end if
       end if
       qsf = epsl/pk
       write (report_unit, fmt="('The safety factor q =      ',f7.4)") qsf
       write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") shat
       if (abs(shat) <= 1.e-5) then
          write (report_unit, fmt="('This is effectively zero; periodic boundary conditions are assumed.')")
       end if
       write (report_unit, fmt="('and epsilon == r/R = ',f7.4)") eps
       write (report_unit, *)
       if (eps > epsilon(0.0)) then
          write (report_unit, fmt="('Trapped particles are included.')")
       else
          write (report_unit, fmt="('Trapped particles are neglected.')")
       end if
       write (report_unit, *)

       if (shift > -epsilon(0.0)) then
          write (report_unit, fmt="('The s-alpha alpha parameter is ',f7.4)") shift
          !CMR 10/11/06: correct sign of dbeta/drho in s-alpha
          write (report_unit, fmt="('corresponding to d beta / d rho = ',f10.4)") -shift/arat/qsf**2
          !CMR 10/11/06: correct sign of dbeta/drho in s-alpha in this check
          if (abs(dbetadrho + shift/arat/qsf**2) > 1.e-2) then
             write (report_unit, *)
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, fmt="('This is inconsistent with beta and the pressure gradient.')")
             write (report_unit, fmt="('################# WARNING #######################')")
          end if
       else
          write (report_unit, *)
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('The s-alpha alpha parameter is less that zero.')")
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
          write (report_unit, fmt="('################# WARNING #######################')")
       end if

    else
       arat = 1.
       write (report_unit, &
            & fmt="('The radius of curvature is infinite.  This is a slab calculation.')")
    end if

    write (report_unit, *)
    select case (model_switch)

    case (model_salpha,model_b2,model_eps)
       if (epsl > 0.) then
          write (report_unit, fmt="('An s-alpha model equilibrium has been selected.')")
          write (report_unit, fmt="('The curvature and grad-B drifts are equal.')")
          write (report_unit, *)
          if (model_switch /= model_eps) then
             write (report_unit, fmt="('For theta0 = 0, each is of the form')")
             write (report_unit, *)
             write (report_unit, fmt="('  epsl*(cos(theta) + (shat*theta-shift*sin(theta))*sin(theta))')")
             write (report_unit, *)
          else
             write (report_unit, fmt="('For theta0 = 0, each is of the form')")
             write (report_unit, *)
             write (report_unit, fmt="('  epsl*(cos(theta) - eps + (shat*theta-shift*sin(theta))*sin(theta))')")
             write (report_unit, *)
          end if
          write (report_unit, fmt="('For finite theta0, there is also a term')")
          write (report_unit, *)
          write (report_unit, fmt="('  -epsl*shat*sin(theta)*theta0')")
          write (report_unit, *)
       end if
       write (report_unit, *)
       write (report_unit, fmt="('For theta0 = 0, |(grad S)**2| is of the form')")
       write (report_unit, *)
       write (report_unit, fmt="('  1.0 + (shat*theta-shift*sin(theta))**2')")
       write (report_unit, *)
       write (report_unit, fmt="('For finite theta0, there is also a term')")
       write (report_unit, *)
       write (report_unit, fmt="('  -shat*(shat*theta - shift*sin(theta))*theta0')")
       write (report_unit, *)
       write (report_unit, fmt="('and finally, the term')")
       write (report_unit, *)
       write (report_unit, fmt="('  shat**2 * theta0**2')")
       write (report_unit, *)
       if (model_switch == model_eps) then
          write (report_unit, *)
          write (report_unit, fmt="(' This model differs from the normal s-alpha model')")
          write (report_unit, fmt="(' only in the curv and grad_B drifts.')")
       end if
       if (model_switch == model_b2) then
          write (report_unit, *)
          write (report_unit, fmt="(' This model differs from the normal s-alpha model')")
          write (report_unit, fmt="(' by an additional factor of 1/B(theta)**2 (not shown above)')")
          write (report_unit, fmt="(' in the curv and grad_B drifts.')")
       end if
    case (model_ccurv)
       write (report_unit, fmt="('Constant curvature is assumed.')")
       write (report_unit, fmt="('The grad-B and curvature drifts are each = ',f10.4)") epsl
       write (report_unit, *)
       write (report_unit, fmt="('For theta0 = 0, |(grad S)**2| is of the form')")
       write (report_unit, *)
       write (report_unit, fmt="('  1.0 + (shat*theta-shift*sin(theta))**2')")
       write (report_unit, *)
       write (report_unit, fmt="('For finite theta0, there is also a term')")
       write (report_unit, *)
       write (report_unit, fmt="('  -shat*shat*theta*theta0')")
       write (report_unit, *)
       write (report_unit, fmt="('and finally, the term')")
       write (report_unit, *)
       write (report_unit, fmt="('  shat**2 * theta0**2')")
       write (report_unit, *)
    case (model_nocurve)
       write (report_unit, fmt="('Zero curvature is assumed.')")
       write (report_unit, *)
       write (report_unit, fmt="('For theta0 = 0, |(grad S)**2| is of the form')")
       write (report_unit, *)
       write (report_unit, fmt="('  1.0 + (shat*theta)**2')")
       write (report_unit, *)
       write (report_unit, fmt="('For finite theta0, there is also a term')")
       write (report_unit, *)
       write (report_unit, fmt="('  -shat*shat*theta*theta0')")
       write (report_unit, *)
       write (report_unit, fmt="('and finally, the term')")
       write (report_unit, *)
       write (report_unit, fmt="('  shat**2 * theta0**2')")
       write (report_unit, *)
    end select
  end subroutine check_theta_grid_salpha