diff --git a/input/equil_a.txt b/input/equil_a.txt index 2685e87..955f442 100644 --- a/input/equil_a.txt +++ b/input/equil_a.txt @@ -1,4 +1,8 @@ - : rr0m,zr0m,rpam ! rhot[m] = min(sqrt((r-rr0m)**2+(z-zr0m)**2), rpam); tor flux phi = pi*b0*rhot**2 - : b0 ! Bphi[T] @ rr0m[m] - : q0,qa,alq ! q = q0 + (qa-q0)*rhot**alq - : nlim ! number of points in first wall (rlim,zlim) polygon +2.96 0.0 1.25 ! R₀,z₀,a, where ρ_p(R, z) = √[(R - R₀)² + (z - z₀)²]/a (m) +6 ! B₀, where B_φ(R) = B₀ R₀/R (T) +3.5 10 2 ! q₀,q₁,α, where q(ρ_p) = q₀ + (q₁-q₀)ρ_p^α +0 ! number of points in the limiter contourn, (R,z) pairs + +! Notes: +! 1. use B₀>0 for clockwise B_φ +! 2. use q₀,q₁<0 for clockwise I_p,B_φ diff --git a/input/profil_a.txt b/input/profil_a.txt index b53543b..04c04dd 100644 --- a/input/profil_a.txt +++ b/input/profil_a.txt @@ -1,3 +1,3 @@ -: dens0,aln1,aln2 ! dens[1e19 m-3] = dens0*(1 - psin**aln1)**aln2 -: te0,dte0,alt1,alt2 ! temp[keV] = (te0-dte0)*(1 - psin**alt1)**alt2 + dte0 -: zeffan +10 0.3 3 ! n₀,a,b, where n(ψ) = n₀(1 - ψ^a)^b (10¹⁹ m⁻³) +14 0 2 8 ! T₀,T₁,a,b, where T(ψ) = (T₀ - T₁)(1 - ψ^a)^b + T₁ (keV) +1.0 ! Z_eff diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index 4708688..65195f9 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -13,9 +13,13 @@ module equilibrium ! Parameters of the analytical equilibrium model type analytic_model - real(wp_) :: q0 ! Safety factor at the magnetic axis - real(wp_) :: q1 ! Safety factor at the edge - real(wp_) :: lq ! Exponent for the q(ρ) power law + real(wp_) :: q0 ! Safety factor at the magnetic axis + real(wp_) :: q1 ! Safety factor at the edge + real(wp_) :: alpha ! Exponent for the q(ρ_p) power law + real(wp_) :: R0 ! R of the magnetic axis (m) + real(wp_) :: z0 ! z of the magnetic axis (m) + real(wp_) :: a ! Minor radius (m) + real(wp_) :: B0 ! Magnetic field at the magnetic axis (T) end type ! Order of the splines @@ -33,13 +37,18 @@ module equilibrium type(analytic_model), save :: model ! More parameters - real(wp_), save :: fpolas - real(wp_), save :: psia, psiant, psinop - real(wp_), save :: phitedge, aminor - real(wp_), save :: btaxis,rmaxis,zmaxis,sgnbphi - real(wp_), save :: btrcen ! used only for jcd_astra def. - real(wp_), save :: rcen ! computed as fpol(a)/btrcen - real(wp_), save :: rmnm, rmxm, zmnm, zmxm + real(wp_), save :: fpolas ! Poloidal current at the edge (r=a), F_a + real(wp_), save :: psia ! Poloidal flux at the edge (r=a), ψ_a + real(wp_), save :: psiant ! Normalised poloidal flux at the true edge + real(wp_), save :: psinop ! Normalised poloidal flux at the true O point + real(wp_), save :: phitedge ! Toroidal flux at the edge + real(wp_), save :: btaxis ! B₀: B_φ = B₀R₀/R + real(wp_), save :: rmaxis, zmaxis ! R,z of the magnetic axis (O point) + real(wp_), save :: sgnbphi ! Sign of B_φ (>0 means CCW) + real(wp_), save :: btrcen ! used only for jcd_astra def. + real(wp_), save :: rcen ! Major radius R₀, computed as fpolas/btrcen + real(wp_), save :: rmnm, rmxm ! R interval of the equilibrium definition + real(wp_), save :: zmnm, zmxm ! z interval of the equilibrium definition real(wp_), save :: zbinf, zbsup private @@ -48,17 +57,17 @@ module equilibrium public equian ! Analytical model public equinum_psi, equinum_fpol ! Interpolated data public fq, bfield, tor_curr, tor_curr_psi ! Accessing local quantities - public frhopol, frhopolv, frhotor ! Converting between poloidal/toroidal flux + public frhotor, frhopol ! Converting between poloidal/toroidal flux public set_equil_spline, set_equil_an ! Initialising internal state public unset_equil_spline ! Deinitialising internal state ! Members exposed for magsurf_data - public kspl, psi_spline, q_spline, points_tgo + public kspl, psi_spline, q_spline, points_tgo, model ! Members exposed to gray_core and more public psia, psiant, psinop - public phitedge, aminor - public btaxis,rmaxis,zmaxis,sgnbphi + public phitedge + public btaxis, rmaxis, zmaxis, sgnbphi public btrcen, rcen public rmnm, rmxm, zmnm, zmxm public zbinf, zbsup @@ -171,7 +180,7 @@ contains ! End of G-EQDSK file close(u) - + ! Build rv,zv,psinr arrays zleft = zmid-0.5_wp_*dz dr = dr/(nr-1) @@ -213,7 +222,6 @@ contains ! local variables integer :: i, u, nlim - real(wp_) :: rr0m, zr0m, rpam, b0 u = get_free_unit(unit) @@ -225,20 +233,9 @@ contains return end if - read(u, *) rr0m, zr0m, rpam - read(u, *) b0 - read(u, *) model%q0, model%q1, model%lq - - if(allocated(data%rv)) deallocate(data%rv) - if(allocated(data%zv)) deallocate(data%zv) - if(allocated(data%fpol)) deallocate(data%fpol) - if(allocated(data%qpsi)) deallocate(data%qpsi) - allocate(data%rv(2), data%zv(1), data%fpol(1), data%qpsi(3)) - - data%rv = [rr0m, rpam] - data%zv = [zr0m] - data%fpol = [b0*rr0m] - data%qpsi = [model%q0, model%q1, model%lq] + read(u, *) model%R0, model%z0, model%a + read(u, *) model%B0 + read(u, *) model%q0, model%q1, model%alpha if(ipass >= 2) then ! get size of boundary and limiter arrays and allocate them @@ -252,7 +249,6 @@ contains read(u,*) (data%rlim(i), data%zlim(i), i = 1, nlim) end if end if - close(u) end subroutine read_equil_an @@ -340,18 +336,19 @@ contains subroutine scale_equil(params, data) - ! Rescale the magnetic field (B) and the plasma current (I) + ! Rescale the magnetic field (B) and the plasma current (I_p) ! and/or force their signs. ! ! Notes: ! 1. signi and signb are ignored on input if equal to 0. - ! They are used to assign the direction of Bphi and Ipla BEFORE scaling. - ! 2. cocos=3 assumed: CCW direction is >0 - ! 3. Bphi and Ipla scaled by the same factor factb to keep q unchanged + ! They are used to assign the direction of B_φ and I_p BEFORE scaling. + ! 2. cocos=3 assumed: positive toroidal direction is CCW from above + ! 3. B_φ and I_p scaled by the same factor factb to keep q unchanged ! 4. factb<0 reverses the directions of Bphi and Ipla use const_and_precisions, only : one use gray_params, only : equilibrium_parameters, equilibrium_data + use gray_params, only : iequil implicit none @@ -359,20 +356,42 @@ contains type(equilibrium_parameters), intent(inout) :: params type(equilibrium_data), intent(inout) :: data - ! local variables - real(wp_) :: last_fpol + ! Notes on cocos=3 + ! 1. With both I_p,B_φ > 0 we have ∂ψ/∂r<0 and ∂Φ/∂r>0. + ! 2. ψ_a = ψ_edge - ψ_axis < 0. + ! 3. q = 1/2π ∂Φ/∂ψ ~ ∂Φ/∂r⋅∂r/∂ψ < 0. + ! 4. In general, sgn(q) = -sgn(I_p)⋅sgn(B_φ). - last_fpol = data%fpol(size(data%fpol)) - - if (params%sgni /=0) & - data%psia = sign(data%psia, real(-params%sgni, wp_)) - if (params%sgnb /=0 .and. params%sgnb * last_fpol < 0) & - data%fpol = -data%fpol + if (iequil < 2) then + ! Analytical model - data%psia = data%psia * params%factb - data%fpol = data%fpol * params%factb - params%sgni = int(sign(one, -data%psia)) - params%sgnb = int(sign(one, last_fpol)) + ! Apply signs + if (params%sgni /= 0) then + model%q0 = sign(model%q0, -params%sgni*one) + model%q1 = sign(model%q1, -params%sgni*one) + end if + if (params%sgnb /= 0) then + model%B0 = sign(model%B0, +params%sgnb*one) + end if + ! Rescale + model%B0 = model%B0 * params%factb + else + ! Numeric data + + ! Apply signs + if (params%sgni /= 0) & + data%psia = sign(data%psia, -params%sgni*one) + if (params%sgnb /= 0) & + data%fpol = sign(data%fpol, +params%sgnb*one) + ! Rescale + data%psia = data%psia * params%factb + data%fpol = data%fpol * params%factb + + ! Compute the signs to be shown in the outputs header when cocos≠0,10. + ! Note: In these cases the values sgni,sgnb from gray.ini are unused. + params%sgni = int(sign(one, -data%psia)) + params%sgnb = int(sign(one, +data%fpol(size(data%fpol)))) + end if end subroutine scale_equil @@ -713,7 +732,7 @@ contains call q_spline%init(psinq, q, size(q)) - ! Toroidal flux φ = 2π ∫q(ψ)dψ + ! Toroidal flux as Φ(ψ) = 2π ∫q(ψ)dψ phit(1)=0 do k=1,q_spline%ndata-1 dx=q_spline%data(k+1)-q_spline%data(k) @@ -726,83 +745,32 @@ contains end subroutine setqphi_num - subroutine set_equil_an(data, n) + subroutine set_equil_an ! Computes the analytical equilibrium data and stores them ! in their respective global variables, see the top of this file. - use const_and_precisions, only : pi, zero, one - use gray_params, only : equilibrium_data + use const_and_precisions, only : pi, one implicit none - ! subroutine arguments - type(equilibrium_data), intent(in) :: data - integer, optional, intent(in) :: n + real(wp_) :: dq - ! local variables - integer, parameter :: nqdef=101 - integer :: i - real(wp_) :: dr, fq0, fq1, qq, res, rn - real(wp_) :: rax, zax, a, bax, qax, q1, qexp - real(wp_), dimension(:), allocatable :: rhotn,rhopn + btaxis = model%B0 + rmaxis = model%R0 + zmaxis = model%z0 + btrcen = model%B0 + rcen = model%R0 + zbinf = model%z0 - model%a + zbsup = model%z0 + model%a + sgnbphi = sign(one, model%B0) - rax = data%rv(1) - zax = data%zv(1) - a = data%rv(2) - bax = data%fpol(1) / rax - qax = data%qpsi(1) - q1 = data%qpsi(2) - qexp = data%qpsi(3) - - btaxis = bax - rmaxis = rax - zmaxis = zax - btrcen = bax - rcen = rax - aminor = a - zbinf = zmaxis-a - zbsup = zmaxis+a - model%q0 = qax - model%q1 = q1 - model%lq = qexp - sgnbphi = sign(one, bax) - - rmxm = rmaxis+aminor - rmnm = rmaxis-aminor + rmxm = model%R0 + model%a + rmnm = model%R0 - model%a zmxm = zbsup zmnm = zbinf - if (present(n)) then - q_spline%ndata = n - else - q_spline%ndata = nqdef - end if - - if (allocated(q_spline%data)) deallocate(q_spline%data) - allocate(q_spline%data(q_spline%ndata)) - allocate(rhotn(q_spline%ndata), rhopn(q_spline%ndata)) - - dr = one/(q_spline%ndata - 1) - rhotn(1) = zero - q_spline%data(1) = zero - res = zero - fq0 = zero - do i = 2, q_spline%ndata - rn = (i - 1)*dr - qq = model%q0 + (model%q1 - model%q0) * rn**model%lq - fq1 = rn / qq - res = res + (fq1 + fq0)/2 * dr - fq0 = fq1 - rhotn(i) = rn - q_spline%data(i) = res - end do - - phitedge = btaxis*aminor**2 ! temporary - psia = res*phitedge - phitedge = pi*phitedge ! final - q_spline%data = q_spline%data/res - rhopn = sqrt(q_spline%data) - - call set_rho_spline(rhopn, rhotn) + phitedge = model%B0 * pi * model%a**2 + dq = (model%q1 - model%q0) / (model%alpha/2 + 1) + psia = 1/(2*pi) * phitedge / (model%q0 + dq) end subroutine set_equil_an @@ -874,120 +842,142 @@ contains end subroutine equinum_fpol - subroutine equian(R, z, psi, fpolv, dfpv, dpsidr, dpsidz, & + subroutine equian(R, z, psin, fpolv, dfpv, dpsidr, dpsidz, & ddpsidrr, ddpsidzz, ddpsidrz) ! Computes all MHD equilibrium parameters from a simple analytical model implicit none ! subroutine arguments - real(wp_), intent(in) :: R,z - real(wp_), intent(out), optional :: psi, fpolv, dfpv, dpsidr, dpsidz, & - ddpsidrr, ddpsidzz, ddpsidrz + real(wp_), intent(in) :: R, z + real(wp_), intent(out), optional :: psin, fpolv, dfpv, dpsidr, dpsidz, & + ddpsidrr, ddpsidzz, ddpsidrz - ! local variables - real(wp_) :: cst, dpsidrp, d2psidrp, dqq, qq, & - rn, rpm, snt, rhop, rhot, btaxqq + ! variables + real(wp_) :: r_p ! poloidal radius + real(wp_) :: rho_p ! poloidal radius normalised - ! simple model for equilibrium: large aspect ratio - ! outside plasma: analytical continuation, not solution Maxwell eqs + ! ρ_p is just the geometrical poloidal radius divided by a + r_p = hypot(R - model%R0, z - model%z0) + rho_p = r_p / model%a - ! Note: rpm is ρ_t in meters, rn is ρ_t normalised - rpm = hypot(R - rmaxis, z - zmaxis) - rn = rpm/aminor + ! ψ_n = ρ_p²(R,z) = [(R-R₀)² + (z-z₀)²]/a² + if (present(psin)) psin = rho_p**2 - snt = 0 - cst = 1 - if (rpm > 0) then - snt = (z - zmaxis)/rpm - cst = (R - rmaxis)/rpm - end if + ! ∂ψ/∂R, ∂ψ/∂z, note that ψ = ψ_n⋅ψ_a + if (present(dpsidr)) dpsidr = 2*(R - model%R0)/model%a**2 * psia + if (present(dpsidz)) dpsidz = 2*(z - model%z0)/model%a**2 * psia - if (present(psi)) then - rhot=rn - if(rn <= 1) then - rhop = frhopol(rhot) - psi = rhop**2 - else - psi = 1 + btaxis / (2*psia*model%q1) * (rpm**2 - aminor**2) - rhop = sqrt(psi) - end if - end if + ! ∂²ψ/∂R², ∂²ψ/∂z², ∂²ψ/∂R∂z + if (present(ddpsidrr)) ddpsidrr = 2/model%a**2 * psia + if (present(ddpsidzz)) ddpsidzz = 2/model%a**2 * psia + if (present(ddpsidrz)) ddpsidrz = 0 - if(rn <= 1) then - qq = model%q0 + (model%q1 - model%q0) * rn**model%lq - btaxqq = btaxis/qq - dpsidrp = btaxqq*rpm - dqq = model%lq * (model%q1 - model%q0) * rn**(model%lq - 1) - d2psidrp=btaxqq*(1 - rn*dqq/qq) - else - btaxqq = btaxis / model%q1 - dpsidrp = btaxqq * rpm - d2psidrp = btaxqq - end if - - if(present(fpolv)) fpolv = btaxis * rmaxis - if(present(dfpv)) dfpv = 0 - - if(present(dpsidr)) dpsidr = dpsidrp*cst - if(present(dpsidz)) dpsidz = dpsidrp*snt - if(present(ddpsidrr)) ddpsidrr = btaxqq*snt**2 + cst**2*d2psidrp - if(present(ddpsidrz)) ddpsidrz = cst * snt * (d2psidrp - btaxqq) - if(present(ddpsidzz)) ddpsidzz = btaxqq * cst**2 + snt**2 * d2psidrp + ! F(ψ) = B₀⋅R₀, a constant + if (present(fpolv)) fpolv = model%B0 * model%R0 + if (present(dfpv)) dfpv = 0 end subroutine equian - function frhopol(rhot) - ! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius - - implicit none - - ! function arguments - real(wp_), intent(in) :: rhot - real(wp_) :: frhopol - - frhopol = rhop_spline%eval(rhot) - end function frhopol - - - function frhopolv(rhot) - ! Vector variant of `frhopol` - use utils, only : locate - - implicit none - - ! function arguments - real(wp_), intent(in) :: rhot(:) - real(wp_) :: frhopolv(size(rhot)) - - ! local variables - integer :: i, i0, j - - i0 = 1 - do j = 1, size(rhot) - call locate(rhop_spline%xdata(i0:), rhop_spline%ndata-i0+1, rhot(j), i) - i = min(max(1,i), rhop_spline%ndata - i0) + i0 - 1 - frhopolv(j) = rhop_spline%raw_eval(i, rhot(j)) - i0 = i - end do - end function frhopolv - - - function frhotor(rhop) + function frhotor(rho_p) ! Converts from poloidal (ρ_p) to toroidal (ρ_t) normalised radius + use gray_params, only : iequil implicit none ! function arguments - real(wp_), intent(in) :: rhop + real(wp_), intent(in) :: rho_p real(wp_) :: frhotor - frhotor = rhot_spline%eval(rhop) + if (iequil < 2) then + ! Analytical model + block + ! The change of variable is obtained by integrating + ! + ! q(ψ) = 1/2π ∂Φ/∂ψ + ! + ! and defining ψ = ψ_a ρ_p², Φ = Φ_a ρ_t². + ! The result is: + ! + ! - ψ_a = 1/2π Φ_a / [q₀ + Δq] + ! + ! - ρ_t = ρ_p √[(q₀ + Δq ρ_p^α)/(q₀ + Δq)] + ! + ! where Δq = (q₁ - q₀)/(α/2 + 1) + real(wp_) :: dq + dq = (model%q1 - model%q0) / (model%alpha/2 + 1) + frhotor = rho_p * sqrt((model%q0 + dq*rho_p**model%alpha) & + / (model%q0 + dq)) + end block + else + ! Numerical data + frhotor = rhot_spline%eval(rho_p) + end if end function frhotor + function frhopol(rho_t) + ! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius + use gray_params, only : iequil + use const_and_precisions, only : comp_eps + + implicit none + + ! function arguments + real(wp_), intent(in) :: rho_t + real(wp_) :: frhopol + + if (iequil < 2) then + ! Analytical model + block + ! In general there is no closed form for ρ_p(ρ_t) in the + ! analytical model, we thus solve numerically the equation + ! ρ_t(ρ_p) = ρ_t₀ for ρ_p. + use minpack, only : hybrj1 + + real(wp_) :: rho_p(1), fvec(1), fjac(1,1), wa(7) + integer :: info + + rho_p = [rho_t] ! first guess, ρ_p ≈ ρ_t + call hybrj1(equation, n=1, x=rho_p, fvec=fvec, fjac=fjac, & + ldfjac=1, tol=comp_eps, info=info, wa=wa, lwa=7) + frhopol = rho_p(1) + end block + else + ! Numerical data + frhopol = rhop_spline%eval(rho_t) + end if + + contains + + subroutine equation(n, x, f, df, ldf, flag) + ! The equation to solve: f(x) = ρ_t(x) - ρ_t₀ = 0 + implicit none + + ! optimal step size + real(wp_), parameter :: e = comp_eps**(1/3.0_wp_) + + ! subroutine arguments + integer, intent(in) :: n, ldf, flag + real(wp_), intent(in) :: x(n) + real(wp_), intent(inout) :: f(n), df(ldf,n) + + if (flag == 1) then + ! return f(x) + f(1) = frhotor(x(1)) - rho_t + else + ! return f'(x), computed numerically + df(1,1) = (frhotor(x(1) + e) - frhotor(x(1) - e)) / (2*e) + end if + end subroutine + + end function frhopol + + function fq(psin) - ! Computes the safety factor q as a function of the - ! normalised poloidal flux ψ + ! Computes the safety factor q as a function of the normalised + ! poloidal flux ψ. + ! + ! Note: this returns the absolute value of q. use gray_params, only : iequil implicit none @@ -996,14 +986,17 @@ contains real(wp_), intent(in) :: psin real(wp_) :: fq - ! local variables - real(wp_) :: rn - if (iequil < 2) then - ! q = q₀ + (q₁ - q₀)ρ^l - rn = frhotor(sqrt(psin)) - fq = model%q0 + (model%q1 - model%q0) * rn**model%lq + ! Analytical model + ! The safety factor is a power law in ρ_p: + ! q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α + block + real(wp_) :: rho_p + rho_p = sqrt(psin) + fq = abs(model%q0 + (model%q1 - model%q0) * rho_p**model%alpha) + end block else + ! Numerical data fq = q_spline%eval(psin) end if end function fq @@ -1024,9 +1017,11 @@ contains real(wp_) :: psin,fpol if (iequil < 2) then + ! Analytical model call equian(rpsim,zpsim,fpolv=bphi,dpsidr=bz,dpsidz=br) if (present(bphi)) bphi=bphi/rpsim else + ! Numerical data call equinum_psi(rpsim,zpsim,psi=bphi,dpsidr=bz,dpsidz=br) if (present(bphi)) then psin=bphi @@ -1034,6 +1029,7 @@ contains bphi=fpol/rpsim end if end if + if (present(br)) br=-br/rpsim if (present(bz)) bz= bz/rpsim end subroutine bfield @@ -1053,14 +1049,16 @@ contains ! local variables real(wp_) :: bzz,dbvcdc13,dbvcdc31 real(wp_) :: dpsidr,ddpsidrr,ddpsidzz - + if(iequil < 2) then + ! Analytical model call equian(rpsim,zpsim,dpsidr=dpsidr, & ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) - else + else + ! Numerical data call equinum_psi(rpsim,zpsim,dpsidr=dpsidr, & ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) - end if + end if bzz= dpsidr/rpsim dbvcdc13=-ddpsidzz/rpsim dbvcdc31= ddpsidrr/rpsim-bzz/rpsim @@ -1105,10 +1103,12 @@ contains character(64) :: msg if (iequil < 2) then - val=frhotor(sqrt(psin)) - r1=rmaxis-val*aminor - r2=rmaxis+val*aminor + ! Analytical model + val = sqrt(psin) + r1 = model%R0 - val*model%a + r2 = model%R0 + val*model%a else + ! Numerical data iopt=1 zc=zmaxis call profil(iopt, psi_spline%knots_x, psi_spline%nknots_x, & @@ -1125,7 +1125,7 @@ contains czc, zeroc, mest, m, ier) r1=zeroc(1) r2=zeroc(2) - end if + end if end subroutine psi_raxis @@ -1187,7 +1187,7 @@ contains case(1) call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz) fvec(1) = dpsidr/psia - fvec(2) = dpsidz/psia + fvec(2) = dpsidz/psia case(2) call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, & ddpsidrz=ddpsidrz) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 5f3468b..e902540 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -2328,7 +2328,7 @@ bb: do do k = 1, q_spline%ndata zk = z(k) if (iequil < 2) then - call equian(rj, zk, psi=psin, fpolv=bphi, dpsidr=bz, dpsidz=br) + call equian(rj, zk, psin=psin, fpolv=bphi, dpsidr=bz, dpsidz=br) else call equinum_psi(rj, zk, psi=psin, dpsidr=bz, dpsidz=br) call equinum_fpol(psin, fpol=bphi) diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index 5509d83..d7fced9 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -447,7 +447,7 @@ contains use gray_params, only : iequil use logger, only : log_warning use dierckx, only : profil, sproota - use equilibrium, only : rmaxis, zmaxis, aminor, frhotor, psi_spline, & + use equilibrium, only : model, frhotor, psi_spline, & kspl, psiant, psinop, points_tgo use limiter, only : rwallm @@ -463,7 +463,7 @@ contains ! local variables integer :: npoints,np,info,ic,ier,iopt,m - real(wp_) :: ra,rb,za,zb,rn,th,zc,val + real(wp_) :: ra,rb,za,zb,th,zc,val real(wp_), dimension(mest) :: zeroc real(wp_), dimension(psi_spline%nknots_x) :: czc @@ -472,13 +472,14 @@ contains th=pi/dble(np) if (iequil<2) then - - rn=frhotor(sqrt(h)) - do ic=1,npoints - zcn(ic)=zmaxis+aminor*rn*sin(th*(ic-1)) - rcn(ic)=rmaxis+aminor*rn*cos(th*(ic-1)) - end do - + block + real(wp_) :: r_p ! poloidal radius + r_p = sqrt(h) * model%a + do ic=1,npoints + rcn(ic) = model%R0 + r_p * cos(th*(ic-1)) + zcn(ic) = model%z0 + r_p * sin(th*(ic-1)) + end do + end block else ra=rup diff --git a/src/main.f90 b/src/main.f90 index c4f9730..95c3de8 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -1,5 +1,5 @@ program main - use const_and_precisions, only : wp_, one, zero + use const_and_precisions, only : wp_, zero use logger, only : INFO, ERROR, set_log_level, log_message use units, only : set_active_units, close_units use utils, only : dirname @@ -209,7 +209,7 @@ contains ! GRAY parameters and data. use equilibrium, only : read_equil_an, read_eqdsk, change_cocos, & set_equil_an, set_equil_spline, scale_equil - use logger, only : log_debug + use logger, only : log_debug implicit none @@ -226,11 +226,6 @@ contains params%raytracing%ipass, & data%equilibrium, err) if (err /= 0) return - - ! Set psia sign to give the correct sign to Iphi - ! (COCOS=3: psia<0 for Iphi>0) - data%equilibrium%psia = sign(one, data%equilibrium%qpsi(2) & - * data%equilibrium%fpol(1)) else ! Numerical equilibrium call log_debug('loading G-EQDK file', & @@ -245,7 +240,7 @@ contains ! Set global variables (for splines) if (params%equilibrium%iequil < 2) then - call set_equil_an(data%equilibrium) + call set_equil_an else call log_debug('computing splines...', mod='main', proc='init_equilibrium') call set_equil_spline(params%equilibrium, data%equilibrium, err) @@ -281,7 +276,7 @@ contains ! temperature, density and plasma effective charge) and initialises ! the respective GRAY data structure. use gray_params, only : profiles_parameters, profiles_data - use equilibrium, only : frhopolv + use equilibrium, only : frhopol use coreprofiles, only : read_profiles_an, read_profiles, & scale_profiles, set_profiles_an, & set_profiles_spline @@ -296,6 +291,9 @@ contains type(profiles_data), intent(out) :: data integer, intent(out) :: err + ! local variables + integer :: i + if (params%iprof == 0) then ! Analytical profiles call log_debug('loading analytical file', & @@ -312,7 +310,7 @@ contains ! Convert psrad to ψ select case (params%irho) case (0) ! psrad is ρ_t = √Φ (toroidal flux) - data%psrad = frhopolv(data%psrad)**2 + data%psrad = [(frhopol(data%psrad(i))**2, i=1,size(data%psrad))] case (1) ! psrad is ρ_p = √ψ (poloidal flux) data%psrad = data%psrad**2 case default ! psrad is already ψ @@ -420,7 +418,10 @@ contains ! Max radius, either due to the plasma extent or equilibrium grid if (params%equilibrium%iequil < 2) then ! analytic equilibrium, use R₀+a - R_max = data%equilibrium%rv(1) + data%equilibrium%rv(2) + block + use equilibrium, only : model + R_max = model%R0 + model%a + end block else ! numeric equilibrium, use max R of the grid R_max = data%equilibrium%rv(size(data%equilibrium%rv))