FIXME : Add documentation
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(in) | :: | report_unit | |||
real, | intent(in) | :: | dbetadrho |
subroutine check_theta_grid_eik(report_unit, dbetadrho)
use theta_grid_params, only: eps
use geometry, only: surf, dp_mult, alpha_input, beta_prime_input, invLp_input
use geometry, only: bishop, irho, eqfile, eqnormfile
use geometry, only: idfit_eq, gen_eq, efit_eq, ppl_eq, local_eq, dfit_eq, chs_eq
implicit none
integer, intent(in) :: report_unit
real, intent(in) :: dbetadrho
real :: beta_prime_new, shat, rhoc
beta_prime_new = eikcoefs_results%dbetadrho
shat = eikcoefs_results%shat
rhoc = eikcoefs_results%rhoc
call checklogic_theta_grid_eik(report_unit)
write (report_unit, *)
if (local_eq) then
write (report_unit, fmt="('The following local equilibrium model has been selected:')")
select case (surf%geoType)
case (0)
write (report_unit, fmt="('Miller')")
case (1)
write (report_unit, fmt="('Global')")
case (2)
write (report_unit, fmt="('GeneralizedEllipticity')")
case (3)
write (report_unit, fmt="('FourierSeries')")
case default
write (report_unit, fmt="('ERROR: invalid analytic geometry specification')")
end select
if (surf%rmaj == 1.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)") surf%rmaj
end if
if (surf%rmaj /= surf%r_geo) then
write (report_unit, *)
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, fmt="('R_geo is not equal to Rmaj.')")
write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
write (report_unit, fmt="('################# WARNING #######################')")
end if
if (irho /= 2) then
write (report_unit, *)
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, fmt="('You have selected irho = ',i2)") irho
write (report_unit, fmt="('For local equilibria, irho=2 is required.')")
write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
write (report_unit, fmt="('################# WARNING #######################')")
end if
write (report_unit, *)
write (report_unit, fmt="('The safety factor q = ',f7.4)") surf%q
eps = rhoc / surf%R_geo
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, *)
write (report_unit, fmt="('B_poloidal is determined by:')")
if (surf%geoType == 1) then
write (report_unit, *)
write (report_unit, fmt="(' minor radius of shaped surface: aSurf = ',f7.4)") surf%aSurf
elseif (surf%geoType == 2 .or. surf%geoType == 3) then
write (report_unit, *)
write (report_unit, fmt="(' first shaping mode number: mMode = ',i2)") surf%mMode
write (report_unit, fmt="(' second shaping mode number: nMode = ',i2)") surf%nMode
end if
write (report_unit, *)
write (report_unit, fmt="(' deltam = ',f7.4)") surf%delm
write (report_unit, fmt="(' deltan = ',f7.4)") surf%deln
if (surf%geoType /= 1) then
write (report_unit, *)
write (report_unit, fmt="(' & gradient: d deltam /d rho = ',f7.4)") surf%delmp
write (report_unit, fmt="(' & gradient: d deltan /d rho = ',f7.4)") surf%delnp
end if
write (report_unit, *)
write (report_unit, fmt="(' thetam = ',f7.4)") surf%thm
write (report_unit, fmt="(' thetan = ',f7.4)") surf%thn
write (report_unit, *)
write (report_unit, fmt="(' shift = ',f7.4)") surf%sHorz
write (report_unit, fmt="(' shiftVert = ',f7.4)") surf%sVert
write (report_unit, *)
write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") shat
write (report_unit, fmt="('This value is set by s_hat_input in the theta_grid_eik_knobs namelist.')")
if (abs(shat) <= 1.e-5) then
write (report_unit, fmt="('This is effectively zero; periodic boundary conditions are assumed.')")
end if
select case (bishop)
case (3)
write (report_unit, fmt="('The normalized inverse pressure gradient scale length = ',f8.4)") invLp_input
case (4)
write (report_unit, fmt="('The beta gradient d beta / d rho = ',f8.4)") beta_prime_input
if (beta_prime_input > epsilon(0.0)) then
write (report_unit, *)
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, fmt="('beta_prime > 0.')")
write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, *)
end if
if (abs(beta_prime_input - dbetadrho) > 1.e-2*max(abs(beta_prime_input),abs(dbetadrho))) then
write (report_unit, *)
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, fmt="('beta_prime_input is not consistent with beta and Lp.')")
write (report_unit, fmt="('beta_prime_input=',f6.4,' dbeta_drho (from species)=',f6.4)") beta_prime_input, dbetadrho
write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, *)
end if
case (5)
write (report_unit, fmt="('The alpha parameter (R beta_prime q**2) = ',f8.4)") alpha_input
! write (*,*) alpha_input, dbetadrho, qinp, Rmaj
if (abs(alpha_input + dbetadrho*surf%q**2*surf%rmaj) > 1.e-2) then
write (report_unit, *)
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, fmt="('alpha is not consistent with beta, q, and Lp.')")
write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, *)
end if
case default
write (report_unit, *)
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, fmt="('You have selected bishop = ',i2)") bishop
write (report_unit, fmt="('For local equilibria, bishop = 4 is recommended.')")
if (bishop == 1) then
write (report_unit, fmt="('For d beta / d rho = 0, bishop = 1 is ok.')")
write (report_unit, fmt="('Otherwise, ')")
end if
write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, *)
end select
end if
if (.not. local_eq) then
if (gen_eq) then
write (report_unit, *)
write (report_unit, fmt="('Equilibrium information obtained from NetCDF file:')")
write (report_unit, fmt="(a)") trim(eqfile)
end if
if (ppl_eq) then
write (report_unit, *)
write (report_unit, fmt="('Equilibrium information obtained from NetCDF file:')")
write (report_unit, fmt="(a)") trim(eqfile)
end if
if (dfit_eq) then
write (report_unit, *)
write (report_unit, fmt="('Dipole equilibrium information obtained from file:')")
write (report_unit, fmt="(a)") trim(eqfile)
write (report_unit, fmt="(a)") trim(eqnormfile)
end if
if (idfit_eq) then
write (report_unit, *)
write (report_unit, fmt="('Dipole equilibrium information obtained from file:')")
write (report_unit, fmt="(a)") trim(eqfile)
end if
if (efit_eq) then
write (report_unit, *)
write (report_unit, fmt="('Equilibrium information obtained from eqdsk:')")
write (report_unit, fmt="(a)") trim(eqfile)
end if
if (chs_eq) then
write (report_unit, *)
write (report_unit, fmt="('Equilibrium information obtained from chease:')")
write (report_unit, fmt="(a)") trim(eqfile)
end if
select case (bishop)
case (1)
write (report_unit, *)
write (report_unit, fmt="('You have set bishop=1, so dp/drho and s_hat will be found from the equilibrium file.')")
write (report_unit, *)
case (3)
write (report_unit, *)
write (report_unit, fmt="('You have set bishop=3.')")
write (report_unit, *)
write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") surf%shat
write (report_unit, fmt="('This value is set by s_hat_input in the theta_grid_eik_knobs namelist.')")
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="('The normalized inverse pressure gradient scale length = ',f8.4)") invLp_input
case (4)
write (report_unit, *)
write (report_unit, fmt="('You have set bishop=4.')")
write (report_unit, *)
write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") surf%shat
write (report_unit, fmt="('This value is set by s_hat_input in the theta_grid_eik_knobs namelist.')")
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="('The beta gradient d beta / d rho = ',f8.4)") beta_prime_input
if (beta_prime_input > epsilon(0.0)) then
write (report_unit, *)
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, fmt="('beta_prime > 0.')")
write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, *)
end if
if (abs(beta_prime_input - dbetadrho) > 1.e-2*abs(dbetadrho)) then
write (report_unit, *)
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, fmt="('beta_prime_input is not consistent with beta and Lp.')")
write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, *)
end if
case (5)
write (report_unit, *)
write (report_unit, fmt="('You have set bishop=5.')")
write (report_unit, *)
write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") surf%shat
write (report_unit, fmt="('This value is set by s_hat_input in the theta_grid_eik_knobs namelist.')")
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="('The alpha parameter (R beta_prime q**2) = ',f8.4)") alpha_input
write (*,*) alpha_input, dbetadrho, surf%q, surf%rmaj
if (abs(alpha_input + dbetadrho*surf%q**2*surf%rmaj) > 1.e-2) then
write (report_unit, *)
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, fmt="('alpha is not consistent with beta, q, and Lp.')")
write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, *)
end if
case (6)
write (report_unit, *)
write (report_unit, fmt="('You have set bishop=6.')")
write (report_unit, *)
write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") surf%shat
write (report_unit, fmt="('This value is set by s_hat_input in the theta_grid_eik_knobs namelist.')")
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="('The value of dp/drho will be found from the equilibrium file.')")
case (7)
write (report_unit, *)
write (report_unit, fmt="('You have set bishop=7.')")
write (report_unit, fmt="('The value of s_hat will be found from the equilibrium file.')")
write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") shat
write (report_unit, fmt="('The value of dp/drho found from the equilibrium file will be multiplied by',f10.4)") dp_mult
write (report_unit, fmt="('to give beta gradient d beta / d rho = ',f8.4)") beta_prime_new
if (abs(beta_prime_new - dbetadrho) > 1.e-2*abs(dbetadrho)) then
write (report_unit, *)
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, fmt="('beta_prime_new is not consistent with beta and Lp.')")
write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, *)
end if
case default
write (report_unit, *)
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, fmt="('You have selected a value for bishop that is not recommended.')")
write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
write (report_unit, fmt="('################# WARNING #######################')")
write (report_unit, *)
end select
end if
end subroutine check_theta_grid_eik