src/gray_tables.f90: fix floating point exception
Skip ψ=0 when computing the flux surfaces contours.
This commit is contained in:
parent
e9e3a4d697
commit
42fcac0726
@ -220,7 +220,7 @@ contains
|
|||||||
if (params%equilibrium%iequil == EQ_VACUUM) return
|
if (params%equilibrium%iequil == EQ_VACUUM) return
|
||||||
|
|
||||||
! Δρ_p for the uniform ρ_p grid
|
! Δρ_p for the uniform ρ_p grid
|
||||||
drho_p = 1.0_wp_ / (n_surf - 1)
|
drho_p = 1.0_wp_ / n_surf
|
||||||
|
|
||||||
! Guess for first contour reconstruction
|
! Guess for first contour reconstruction
|
||||||
R_hi = equil%axis(1)
|
R_hi = equil%axis(1)
|
||||||
@ -228,7 +228,7 @@ contains
|
|||||||
z_hi = equil%axis(2) + (equil%z_boundary(2) - equil%axis(2))/10
|
z_hi = equil%axis(2) + (equil%z_boundary(2) - equil%axis(2))/10
|
||||||
z_lo = equil%axis(2) - (equil%axis(2) - equil%z_boundary(1))/10
|
z_lo = equil%axis(2) - (equil%axis(2) - equil%z_boundary(1))/10
|
||||||
|
|
||||||
do i = 0, n_surf - 1
|
do i = 1, n_surf - 1
|
||||||
rho_p = i * drho_p
|
rho_p = i * drho_p
|
||||||
psi_n = rho_p**2
|
psi_n = rho_p**2
|
||||||
call equil%flux_contour(psi_n, R, z, R_hi, z_hi, R_lo, z_lo)
|
call equil%flux_contour(psi_n, R, z, R_hi, z_hi, R_lo, z_lo)
|
||||||
|
Loading…
Reference in New Issue
Block a user