istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

replace equilibrium module with an object

Browse Source 
Similarly to eb648039 this change replaces the `equilibrium` module with
a new `gray_equil` module providing the same functionality without using
global variables.

  - `read_eqdsk`, `read_equil_an` are replaced by a single `load_equil`
    routine that handles all equilibrium kind (analytical, numerical,
    and vacuum).

  - `scale_equil` is merged into `load_equil`, which besides reading
    the equilibrium from file peforms the rescaling and interpolation based
    on the `gray_parameters` settings and the equilibrium kind.

    To operate on G-EQDSK data specifically, the `change_cocors` and
    `scale_eqdsk` are still available. The numeric equilibrium must then
    be initialised manually by calling equil%init().

  - `set_equil_spline`, `set_equil_an`, `unset_equil_spline`
     are completely removed as the module no longer has any internal state.

  - `fq` is replaced by `equil%safety`; `bfield` by `equil%b_field`;
    `frhotor`, `frhopol` by `equil%pol2tor` and `equil%pol2tor`;
    and the remaining subroutines by other methods of `abstract_equil`
    retaining the old name.

  - the `contours_psi` subroutine is replaced by `equil%flux_contour`,
    with a slightly changed invocation but same functionality.

  - the `gray_data` type is no longer required ans has been removed: all
    the core subroutines now access the input data only though either
    `abstract_equil`, `abstract_plasma` or the `limiter` contour.
This commit is contained in:
Michele Guerini Rocco 2024-08-29 17:16:33 +02:00 committed by rnhmjoj
parent ae80fb4945
commit 166086d369
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
12 changed files with 1984 additions and 1848 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1343
src/equilibrium.f90
View File

File diff suppressed because it is too large Load Diff

194
src/gray_core.f90
View File

@ -7,11 +7,12 @@ module gray_core
contains contains
subroutine gray_main(params, data, plasma, results, error, rhout) subroutine gray_main(params, equil, plasma, limiter, results, error, rhout)
use const_and_precisions, only : zero, one, comp_tiny use const_and_precisions, only : zero, one, comp_tiny
use polarization, only : ellipse_to_field use polarization, only : ellipse_to_field
use types, only : table, wrap use types, only : contour, table, wrap
use gray_params, only : gray_parameters, gray_data, gray_results, EQ_VACUUM use gray_params, only : gray_parameters, gray_results, EQ_VACUUM
use gray_equil, only : abstract_equil
use gray_plasma, only : abstract_plasma use gray_plasma, only : abstract_plasma
use gray_tables, only : init_tables, store_ec_profiles, store_ray_data, & use gray_tables, only : init_tables, store_ec_profiles, store_ray_data, &
store_beam_shape, find_flux_surfaces, & store_beam_shape, find_flux_surfaces, &
@ -29,8 +30,9 @@ contains
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(inout) :: params type(gray_parameters), intent(inout) :: params
type(gray_data), intent(in) :: data class(abstract_equil), intent(in) :: equil
class(abstract_plasma), intent(in) :: plasma class(abstract_plasma), intent(in) :: plasma
type(contour), intent(in) :: limiter
type(gray_results), intent(out) :: results type(gray_results), intent(out) :: results
! Predefined grid for the output profiles (optional) ! Predefined grid for the output profiles (optional)
@ -115,10 +117,10 @@ contains
if (params%equilibrium%iequil /= EQ_VACUUM) then if (params%equilibrium%iequil /= EQ_VACUUM) then
! Initialise the magsurf_data module ! Initialise the magsurf_data module
call compute_flux_averages(params, results%tables) call compute_flux_averages(params, equil, results%tables)
! Initialise the output profiles ! Initialise the output profiles
call pec_init(params%output, rhout) call pec_init(params%output, equil, rhout)
end if end if
! Allocate memory for the results... ! Allocate memory for the results...
@ -133,14 +135,14 @@ contains
! ========= set environment END ========= ! ========= set environment END =========
! Pre-determinted tables ! Pre-determinted tables
results%tables%kinetic_profiles = kinetic_profiles(params, plasma) results%tables%kinetic_profiles = kinetic_profiles(params, equil, plasma)
results%tables%ec_resonance = ec_resonance(params, bres) results%tables%ec_resonance = ec_resonance(params, equil, bres)
results%tables%inputs_maps = inputs_maps(params, plasma, bres, xgcn) results%tables%inputs_maps = inputs_maps(params, equil, plasma, bres, xgcn)
! print ψ surface for q=3/2 and q=2/1 ! print ψ surface for q=3/2 and q=2/1
call find_flux_surfaces( & call find_flux_surfaces(qvals=[1.5_wp_, 2.0_wp_], &
qvals=[1.5_wp_, 2.0_wp_], params=params, & params=params, equil=equil, &
tbl=results%tables%flux_surfaces) tbl=results%tables%flux_surfaces)
! print initial position ! print initial position
write (msg, '("initial position:",3(x,g0.3))') params%antenna%pos write (msg, '("initial position:",3(x,g0.3))') params%antenna%pos
@ -218,14 +220,14 @@ contains
lgcpl1 = zero lgcpl1 = zero
p0ray = p0jk ! * initial beam power p0ray = p0jk ! * initial beam power
call ic_gb(params, anv0, ak0, yw, ypw, stv, xc, du1, & call ic_gb(params, equil, anv0, ak0, yw, ypw, stv, xc, &
gri, ggri, index_rt, results%tables) ! * initial conditions du1, gri, ggri, index_rt, results%tables) ! * initial conditions
do jk=1,params%raytracing%nray do jk=1,params%raytracing%nray
zzm = yw(3,jk)*0.01_wp_ zzm = yw(3,jk)*0.01_wp_
rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_ rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_
if (data%equilibrium%limiter%contains(rrm, zzm)) then ! * start propagation in/outside vessel? if (limiter%contains(rrm, zzm)) then ! * start propagation in/outside vessel?
iow(jk) = 1 ! + inside iow(jk) = 1 ! + inside
else else
iow(jk) = 0 ! + outside iow(jk) = 0 ! + outside
@ -258,7 +260,7 @@ contains
do jk=1,params%raytracing%nray do jk=1,params%raytracing%nray
if(iwait(jk)) cycle ! jk ray is waiting for next pass if(iwait(jk)) cycle ! jk ray is waiting for next pass
stv(jk) = stv(jk) + params%raytracing%dst ! current ray step stv(jk) = stv(jk) + params%raytracing%dst ! current ray step
call rkstep(params, plasma, sox, bres, xgcn, yw(:,jk), & call rkstep(params, equil, plasma, sox, bres, xgcn, yw(:,jk), &
ypw(:,jk), gri(:,jk), ggri(:,:,jk), igrad_b) ypw(:,jk), gri(:,jk), ggri(:,:,jk), igrad_b)
end do end do
! update position and grad ! update position and grad
@ -274,10 +276,11 @@ contains
! compute derivatives with updated gradient and local plasma values ! compute derivatives with updated gradient and local plasma values
xv = yw(1:3,jk) xv = yw(1:3,jk)
anv = yw(4:6,jk) anv = yw(4:6,jk)
call ywppla_upd(params, plasma, xv, anv, gri(:,jk), ggri(:,:,jk), & call ywppla_upd(params, equil, plasma, &
sox, bres, xgcn,ypw(:,jk), psinv, dens, btot, bv, & xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres, &
xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, & xgcn,ypw(:,jk), psinv, dens, btot, bv, xg, yg, &
derdnm, ierrn, igrad_b) derxg, anpl, anpr, ddr, ddi, dersdst, derdnm, &
ierrn, igrad_b)
! update global error code and print message ! update global error code and print message
if(ierrn/=0) then if(ierrn/=0) then
error = ior(error,ierrn) error = ior(error,ierrn)
@ -289,7 +292,7 @@ contains
rrm = sqrt(xv(1)*xv(1)+xv(2)*xv(2))*0.01_wp_ rrm = sqrt(xv(1)*xv(1)+xv(2)*xv(2))*0.01_wp_
ins_pl = (psinv>=zero .and. psinv<params%profiles%psnbnd) ! in/out plasma? ins_pl = (psinv>=zero .and. psinv<params%profiles%psnbnd) ! in/out plasma?
ins_wl = data%equilibrium%limiter%contains(rrm, zzm) ! in/out vessel? ins_wl = limiter%contains(rrm, zzm) ! in/out vessel?
ent_pl = (mod(iop(jk),2) == 0 .and. ins_pl) ! enter plasma ent_pl = (mod(iop(jk),2) == 0 .and. ins_pl) ! enter plasma
ext_pl = (mod(iop(jk),2) == 1 .and. .not.ins_pl) ! exit plasma ext_pl = (mod(iop(jk),2) == 1 .and. .not.ins_pl) ! exit plasma
ent_wl = (mod(iow(jk),2) == 0 .and. ins_wl) ! enter vessel ent_wl = (mod(iow(jk),2) == 0 .and. ins_wl) ! enter vessel
@ -300,8 +303,9 @@ contains
call log_debug(msg, mod='gray_core', proc='gray_main') call log_debug(msg, mod='gray_core', proc='gray_main')
call ellipse_to_field(psipv(parent_index_rt), chipv(parent_index_rt), & ! compute polarisation and couplings call ellipse_to_field(psipv(parent_index_rt), chipv(parent_index_rt), & ! compute polarisation and couplings
ext, eyt) ext, eyt)
call plasma_in(jk, xv, anv, bres, sox, cpl, psipol, chipol, iop, ext, eyt, & call plasma_in(jk, equil, xv, anv, bres, sox, cpl, &
perfect=.not. params%raytracing%ipol & psipol, chipol, iop, ext, eyt, &
perfect=.not. params%raytracing%ipol &
.and. params%antenna%iox == iox & .and. params%antenna%iox == iox &
.and. iop(jk) == 0 .and. ip==1) .and. iop(jk) == 0 .and. ip==1)
@ -361,22 +365,23 @@ contains
else if(ext_pl) then ! ray exits plasma else if(ext_pl) then ! ray exits plasma
write (msg, '(" ray ", g0, " left plasma")') jk write (msg, '(" ray ", g0, " left plasma")') jk
call log_debug(msg, mod='gray_core', proc='gray_main') call log_debug(msg, mod='gray_core', proc='gray_main')
call plasma_out(jk,xv,anv,bres,sox,iop,ext,eyt) call plasma_out(jk, equil, xv, anv, bres, sox, iop, ext, eyt)
end if end if
if(ent_wl) then ! ray enters vessel if(ent_wl) then ! ray enters vessel
iow(jk)=iow(jk)+1 ! * out->in iow(jk)=iow(jk)+1 ! * out->in
else if(ext_wl) then ! ray exit vessel else if(ext_wl) then ! ray exit vessel
call wall_out(jk, data%equilibrium%limiter, ins_pl, xv, anv, & call wall_out(jk, equil, limiter, ins_pl, xv, anv, &
params%raytracing%dst, bres, sox, psipol, chipol, & params%raytracing%dst, bres, sox, psipol, chipol, &
iow, iop, ext, eyt) iow, iop, ext, eyt)
yw(:,jk) = [xv, anv] ! * updated coordinates (reflected) yw(:,jk) = [xv, anv] ! * updated coordinates (reflected)
igrad_b = 0 ! * switch to ray-tracing igrad_b = 0 ! * switch to ray-tracing
call ywppla_upd(params, plasma, xv, anv, gri(:,jk), ggri(:,:,jk), & call ywppla_upd(params, equil, plasma, &
sox, bres, xgcn, ypw(:,jk), psinv, dens, btot, & xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres, &
bv, xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, & xgcn, ypw(:,jk), psinv, dens, btot, bv, xg, &
yg, derxg, anpl, anpr, ddr, ddi, dersdst, &
derdnm, ierrn, igrad_b) ! * update derivatives after reflection derdnm, ierrn, igrad_b) ! * update derivatives after reflection
if(ierrn/=0) then ! * update global error code and print message if(ierrn/=0) then ! * update global error code and print message
error = ior(error,ierrn) error = ior(error,ierrn)
@ -396,8 +401,8 @@ contains
yypnext(:,jk,index_rt) = ypw(:,jk) ! for next pass = reflection point yypnext(:,jk,index_rt) = ypw(:,jk) ! for next pass = reflection point
stnext(jk,index_rt) = stv(jk) ! . starting step for next pass = step after reflection stnext(jk,index_rt) = stv(jk) ! . starting step for next pass = step after reflection
call plasma_in(jk, xv, anv, bres, sox, cpl, psipol, chipol, & ! . ray re-enters plasma after reflection call plasma_in(jk, equil, xv, anv, bres, sox, cpl, & ! . ray re-enters plasma after reflection
iop, ext, eyt, perfect=.false.) psipol, chipol, iop, ext, eyt, perfect=.false.)
if(cpl(1) < etaucr) then ! . low coupled power for O-mode? => de-activate derived rays if(cpl(1) < etaucr) then ! . low coupled power for O-mode? => de-activate derived rays
call turnoffray(jk,ip+1,npass,2*ib-1,iroff) call turnoffray(jk,ip+1,npass,2*ib-1,iroff)
@ -434,10 +439,10 @@ contains
tekev = plasma%temp(psinv) tekev = plasma%temp(psinv)
block block
complex(wp_) :: Npr_warm complex(wp_) :: Npr_warm
call alpha_effj(params%ecrh_cd, plasma, psinv, xg, yg, dens, & call alpha_effj(params%ecrh_cd, equil, plasma, &
tekev, ak0, bres, derdnm, anpl, anpr, sox, & psinv, xg, yg, dens, tekev, ak0, bres, &
Npr_warm, alpha, didp, nharm, nhf, iokhawa, & derdnm, anpl, anpr, sox, Npr_warm, alpha, &
ierrhcd) didp, nharm, nhf, iokhawa, ierrhcd)
anprre = Npr_warm%re anprre = Npr_warm%re
anprim = Npr_warm%im anprim = Npr_warm%im
if(ierrhcd /= 0) then if(ierrhcd /= 0) then
@ -481,7 +486,7 @@ contains
ccci(jk,i:params%raytracing%nstep) = ccci(jk,i-1) ccci(jk,i:params%raytracing%nstep) = ccci(jk,i-1)
psjki(jk,i:params%raytracing%nstep) = psjki(jk,i-1) psjki(jk,i:params%raytracing%nstep) = psjki(jk,i-1)
else else
call store_ray_data(params, results%tables, & call store_ray_data(params, equil, results%tables, &
i, jk, stv(jk), p0ray(jk), xv, psinv, btot, bv, ak0, & i, jk, stv(jk), p0ray(jk), xv, psinv, btot, bv, ak0, &
anpl, anpr, anv, anprim, dens, tekev, alpha, tau, dids, & anpl, anpr, anv, anprim, dens, tekev, alpha, tau, dids, &
nharm, nhf, iokhawa, index_rt, ddr, ddi, xg, yg, derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes) nharm, nhf, iokhawa, index_rt, ddr, ddi, xg, yg, derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes)
@ -583,9 +588,9 @@ contains
results%tables%ec_profiles, rhop_tab, rhot_tab, jphi_beam, & results%tables%ec_profiles, rhop_tab, rhot_tab, jphi_beam, &
jcd_beam, dpdv_beam, currins_beam, pins_beam, ip) jcd_beam, dpdv_beam, currins_beam, pins_beam, ip)
call postproc_profiles(pabs_beam,icd_beam,rhot_tab,dpdv_beam, & call postproc_profiles(equil, pabs_beam, icd_beam, rhot_tab, dpdv_beam, &
jphi_beam, rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip, & jphi_beam, rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip, &
rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) ! *compute profiles width for current beam rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx) ! *compute profiles width for current beam
if (results%tables%summary%active) then if (results%tables%summary%active) then
call results%tables%summary%append([ & call results%tables%summary%append([ &
@ -665,18 +670,20 @@ contains
end subroutine vectinit end subroutine vectinit
subroutine ic_gb(params, anv0c, ak0, ywrk0, ypwrk0, & subroutine ic_gb(params, equil, anv0c, ak0, ywrk0, ypwrk0, &
stv, xc0, du10, gri, ggri, index_rt, & stv, xc0, du10, gri, ggri, index_rt, &
tables) tables)
! beam tracing initial conditions igrad=1 ! beam tracing initial conditions igrad=1
! !!!!!! check ray tracing initial conditions igrad=0 !!!!!! ! !!!!!! check ray tracing initial conditions igrad=0 !!!!!!
use const_and_precisions, only : zero,one,pi,half,two,degree,ui=>im use const_and_precisions, only : zero,one,pi,half,two,degree,ui=>im
use gray_params, only : gray_parameters, gray_tables use gray_params, only : gray_parameters, gray_tables
use gray_equil, only : abstract_equil
use types, only : table use types, only : table
use gray_tables, only : store_ray_data use gray_tables, only : store_ray_data
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
real(wp_), dimension(3), intent(in) :: anv0c real(wp_), dimension(3), intent(in) :: anv0c
real(wp_), intent(in) :: ak0 real(wp_), intent(in) :: ak0
real(wp_), dimension(6, params%raytracing%nray), intent(out) :: ywrk0, ypwrk0 real(wp_), dimension(6, params%raytracing%nray), intent(out) :: ywrk0, ypwrk0
@ -984,7 +991,7 @@ contains
! save step "zero" data ! save step "zero" data
if (present(tables)) & if (present(tables)) &
call store_ray_data(params, tables, & call store_ray_data(params, equil, tables, &
i=0, jk=jk, s=stv(jk), P0=one, pos=xc0(:,k,j), & i=0, jk=jk, s=stv(jk), P0=one, pos=xc0(:,k,j), &
psi_n=-one, B=zero, b_n=[zero,zero,zero], k0=ak0, & psi_n=-one, B=zero, b_n=[zero,zero,zero], k0=ak0, &
N_pl=zero, N_pr=zero, N=ywrk0(:, jk), N_pr_im=zero, & N_pl=zero, N_pr=zero, N=ywrk0(:, jk), N_pr_im=zero, &
@ -998,21 +1005,23 @@ contains
subroutine rkstep(params, plasma, sox, bres, xgcn, y, yp, dgr, ddgr, igrad) subroutine rkstep(params, equil, plasma, &
sox, bres, xgcn, y, yp, dgr, ddgr, igrad)
! Runge-Kutta integrator ! Runge-Kutta integrator
use gray_params, only : gray_parameters use gray_params, only : gray_parameters
use gray_equil, only : abstract_equil
use gray_plasma, only : abstract_plasma use gray_plasma, only : abstract_plasma
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_plasma), intent(in) :: plasma class(abstract_equil), intent(in) :: equil
class(abstract_plasma), intent(in) :: plasma
real(wp_), intent(in) :: bres, xgcn real(wp_), intent(in) :: bres, xgcn
real(wp_), intent(inout) :: y(6) real(wp_), intent(inout) :: y(6)
real(wp_), intent(in) :: yp(6) real(wp_), intent(in) :: yp(6)
real(wp_), intent(in) :: dgr(3) real(wp_), intent(in) :: dgr(3)
real(wp_), intent(in) :: ddgr(3,3) real(wp_), intent(in) :: ddgr(3,3)
integer, intent(in) :: igrad,sox integer, intent(in) :: igrad,sox
! local variables ! local variables
real(wp_), dimension(6) :: k1, k2, k3, k4 real(wp_), dimension(6) :: k1, k2, k3, k4
@ -1030,19 +1039,22 @@ contains
function f(y) function f(y)
real(wp_), intent(in) :: y(6) real(wp_), intent(in) :: y(6)
real(wp_) :: f(6) real(wp_) :: f(6)
call rhs(params, plasma, sox, bres, xgcn, y, dgr, ddgr, f, igrad) call rhs(params, equil, plasma, sox, bres, xgcn, y, dgr, ddgr, f, igrad)
end function end function
end subroutine rkstep end subroutine rkstep
subroutine rhs(params, plasma, sox, bres, xgcn, y, dgr, ddgr, dery, igrad) subroutine rhs(params, equil, plasma, &
sox, bres, xgcn, y, dgr, ddgr, dery, igrad)
! Compute right-hand side terms of the ray equations (dery) ! Compute right-hand side terms of the ray equations (dery)
! used in R-K integrator ! used in R-K integrator
use gray_params, only : gray_parameters use gray_params, only : gray_parameters
use gray_equil, only : abstract_equil
use gray_plasma, only : abstract_plasma use gray_plasma, only : abstract_plasma
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
class(abstract_plasma), intent(in) :: plasma class(abstract_plasma), intent(in) :: plasma
real(wp_), intent(in) :: y(6) real(wp_), intent(in) :: y(6)
real(wp_), intent(in) :: bres, xgcn real(wp_), intent(in) :: bres, xgcn
@ -1057,7 +1069,7 @@ contains
real(wp_), dimension(3,3) :: derbv real(wp_), dimension(3,3) :: derbv
xv = y(1:3) xv = y(1:3)
call plas_deriv(params, plasma, xv, bres, xgcn, dens, btot, & call plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, &
bv, derbv, xg, yg, derxg, deryg) bv, derbv, xg, yg, derxg, deryg)
anv = y(4:6) anv = y(4:6)
call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, & call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, &
@ -1065,18 +1077,21 @@ contains
end subroutine rhs end subroutine rhs
subroutine ywppla_upd(params, plasma, xv, anv, dgr, ddgr, sox, bres, xgcn, dery, & subroutine ywppla_upd(params, equil, plasma, &
psinv, dens, btot, bv, xg, yg, derxg, anpl, anpr, & xv, anv, dgr, ddgr, sox, bres, xgcn, dery, &
ddr, ddi, dersdst, derdnm, error, igrad) psinv, dens, btot, bv, xg, yg, derxg, anpl, &
anpr, ddr, ddi, dersdst, derdnm, error, igrad)
! Compute right-hand side terms of the ray equations (dery) ! Compute right-hand side terms of the ray equations (dery)
! used after full R-K step and grad(S_I) update ! used after full R-K step and grad(S_I) update
use gray_errors, only : raise_error, large_npl use gray_errors, only : raise_error, large_npl
use gray_params, only : gray_parameters use gray_params, only : gray_parameters
use gray_equil, only : abstract_equil
use gray_plasma, only : abstract_plasma use gray_plasma, only : abstract_plasma
! subroutine rguments ! subroutine rguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_plasma), intent(in) :: plasma class(abstract_equil), intent(in) :: equil
class(abstract_plasma), intent(in) :: plasma
real(wp_), intent(in) :: xv(3), anv(3) real(wp_), intent(in) :: xv(3), anv(3)
real(wp_), intent(in) :: dgr(3) real(wp_), intent(in) :: dgr(3)
@ -1096,9 +1111,10 @@ contains
real(wp_), dimension(3,3) :: derbv real(wp_), dimension(3,3) :: derbv
real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_ real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_
call plas_deriv(params, plasma, xv,bres,xgcn,dens,btot,bv,derbv,xg,yg,derxg,deryg,psinv) call plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, &
call disp_deriv(params, anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad, & bv, derbv, xg, yg, derxg, deryg, psinv)
dery,anpl,anpr,ddr,ddi,dersdst,derdnm) call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, bv, derbv, dgr, &
ddgr, igrad, dery, anpl, anpr, ddr, ddi, dersdst, derdnm)
if (abs(anpl) > anplth2) then if (abs(anpl) > anplth2) then
error = raise_error(error, large_npl, 1) error = raise_error(error, large_npl, 1)
@ -1315,15 +1331,14 @@ contains
end subroutine solg3 end subroutine solg3
subroutine plas_deriv(params, plasma, xv, bres, xgcn, dens, btot, bv, & subroutine plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, &
derbv, xg, yg, derxg, deryg, psinv_opt) bv, derbv, xg, yg, derxg, deryg, psinv_opt)
use const_and_precisions, only : zero, cm use const_and_precisions, only : zero, cm
use gray_params, only : gray_parameters, EQ_VACUUM use gray_equil, only : abstract_equil, vacuum
use gray_plasma, only : abstract_plasma use gray_plasma, only : abstract_plasma
use equilibrium, only : psia, pol_flux, pol_curr, sgnbphi
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(in) :: params class(abstract_equil), intent(in) :: equil
class(abstract_plasma), intent(in) :: plasma class(abstract_plasma), intent(in) :: plasma
real(wp_), dimension(3), intent(in) :: xv real(wp_), dimension(3), intent(in) :: xv
real(wp_), intent(in) :: xgcn, bres real(wp_), intent(in) :: xgcn, bres
@ -1352,11 +1367,12 @@ contains
bv = zero bv = zero
derbv = zero derbv = zero
if (params%equilibrium%iequil == EQ_VACUUM) then select type (equil)
! copy optional output type is (vacuum)
if (present(psinv_opt)) psinv_opt = psinv ! copy optional output
return if (present(psinv_opt)) psinv_opt = psinv
end if return
end select
dbtot = zero dbtot = zero
dbv = zero dbv = zero
@ -1374,15 +1390,16 @@ contains
csphi = xx/rr csphi = xx/rr
snphi = yy/rr snphi = yy/rr
bv(1) = -snphi*sgnbphi bv(1) = -snphi*equil%sgn_bphi
bv(2) = csphi*sgnbphi bv(2) = csphi*equil%sgn_bphi
! convert from cm to meters ! convert from cm to meters
zzm = 1.0e-2_wp_*zz zzm = 1.0e-2_wp_*zz
rrm = 1.0e-2_wp_*rr rrm = 1.0e-2_wp_*rr
call pol_flux(rrm, zzm, psinv, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz) call equil%pol_flux(rrm, zzm, psinv, dpsidr, dpsidz, &
call pol_curr(psinv, fpolv, dfpv) ddpsidrr, ddpsidzz, ddpsidrz)
call equil%pol_curr(psinv, fpolv, dfpv)
! copy optional output ! copy optional output
if (present(psinv_opt)) psinv_opt = psinv if (present(psinv_opt)) psinv_opt = psinv
@ -1402,8 +1419,8 @@ contains
! B = f(psi)/R e_phi+ grad psi x e_phi/R ! B = f(psi)/R e_phi+ grad psi x e_phi/R
bphi = fpolv/rrm bphi = fpolv/rrm
brr = -dpsidz*psia/rrm brr = -dpsidz*equil%psi_a/rrm
bzz = +dpsidr*psia/rrm bzz = +dpsidr*equil%psi_a/rrm
! bvc(i) = B_i in cylindrical coordinates ! bvc(i) = B_i in cylindrical coordinates
bvc = [brr, bphi, bzz] bvc = [brr, bphi, bzz]
@ -1414,12 +1431,12 @@ contains
bv(3)=bvc(3) bv(3)=bvc(3)
! dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv) ! dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv)
dbvcdc(1,1) = -ddpsidrz*psia/rrm - brr/rrm dbvcdc(1,1) = -ddpsidrz*equil%psi_a/rrm - brr/rrm
dbvcdc(1,3) = -ddpsidzz*psia/rrm dbvcdc(1,3) = -ddpsidzz*equil%psi_a/rrm
dbvcdc(2,1) = dfpv*dpsidr/rrm - bphi/rrm dbvcdc(2,1) = dfpv*dpsidr/rrm - bphi/rrm
dbvcdc(2,3) = dfpv*dpsidz/rrm dbvcdc(2,3) = dfpv*dpsidz/rrm
dbvcdc(3,1) = +ddpsidrr*psia/rrm - bzz/rrm dbvcdc(3,1) = +ddpsidrr*equil%psi_a/rrm - bzz/rrm
dbvcdc(3,3) = +ddpsidrz*psia/rrm dbvcdc(3,3) = +ddpsidrz*equil%psi_a/rrm
! dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv) ! dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv)
dbvdc(1,1) = dbvcdc(1,1)*csphi - dbvcdc(2,1)*snphi dbvdc(1,1) = dbvcdc(1,1)*csphi - dbvcdc(2,1)*snphi
@ -1746,15 +1763,15 @@ contains
subroutine alpha_effj(params, plasma, psinv, X, Y, density, temperature, & subroutine alpha_effj(params, equil, plasma, psinv, X, Y, density, &
k0, Bres, derdnm, Npl, Npr_cold, sox, Npr, & temperature, k0, Bres, derdnm, Npl, Npr_cold, &
alpha, dIdp, nhmin, nhmax, iokhawa, error) sox, Npr, alpha, dIdp, nhmin, nhmax, iokhawa, error)
! Computes the absorption coefficient and effective current ! Computes the absorption coefficient and effective current
use const_and_precisions, only : pi, mc2=>mc2_ use const_and_precisions, only : pi, mc2=>mc2_
use gray_params, only : ecrh_cd_parameters use gray_params, only : ecrh_cd_parameters
use gray_equil, only : abstract_equil
use gray_plasma, only : abstract_plasma use gray_plasma, only : abstract_plasma
use equilibrium, only : sgnbphi
use dispersion, only : harmnumber, warmdisp use dispersion, only : harmnumber, warmdisp
use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl
use gray_errors, only : negative_absorption, raise_error use gray_errors, only : negative_absorption, raise_error
@ -1767,8 +1784,9 @@ contains
! ECRH & CD parameters ! ECRH & CD parameters
type(ecrh_cd_parameters), intent(in) :: params type(ecrh_cd_parameters), intent(in) :: params
! plasma object ! MHD equilibrium, plasma object
class(abstract_plasma), intent(in) :: plasma class(abstract_equil), intent(in) :: equil
class(abstract_plasma), intent(in) :: plasma
! poloidal flux ψ ! poloidal flux ψ
real(wp_), intent(in) :: psinv real(wp_), intent(in) :: psinv
! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω
@ -1899,7 +1917,7 @@ contains
! current drive efficiency <j/p> [A⋅m/W] ! current drive efficiency <j/p> [A⋅m/W]
effjcdav = rbavi*effjcd effjcdav = rbavi*effjcd
dIdp = sgnbphi * effjcdav / (2*pi*rrii) dIdp = equil%sgn_bphi * effjcdav / (2*pi*rrii)
end subroutine alpha_effj end subroutine alpha_effj

1620
src/gray_equil.f90 Normal file
View File

File diff suppressed because it is too large Load Diff

67
src/gray_jetto1beam.f90
View File

@ -1,10 +1,11 @@
subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
nbnd, rbnd, zbnd, nrho, psrad, fpol, te, dne, zeff, qpsi, ibeam, & nbnd, rbnd, zbnd, nrho, psrad, fpol, te, dne, zeff, qpsi, ibeam, &
p0mw, alphain, betain, dpdv, jcd, pabs, icd, err) p0mw, alphain, betain, dpdv, jcd, pabs, icd, err)
use const_and_precisions, only: wp_ use const_and_precisions, only : wp_
use gray_params, only: gray_parameters, gray_data, gray_results use gray_params, only : gray_parameters, gray_results
use gray_core, only: gray_main use gray_core, only : gray_main
use gray_plasma, only: numeric_plasma use gray_equil, only : numeric_equil, contour
use gray_plasma, only : numeric_plasma
implicit none implicit none
@ -21,11 +22,12 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
integer, intent(out) :: err integer, intent(out) :: err
! local variables ! local variables
type(gray_parameters) :: params type(gray_parameters), save :: params
type(gray_data) :: data type(numeric_equil) :: equil
type(numeric_plasma) :: plasma type(numeric_plasma) :: plasma
type(gray_results) :: res type(contour), save :: limiter
logical, save :: firstcall = .true. type(gray_results) :: res
logical, save :: firstcall = .true.
! Initialisation tasks for the first call only ! Initialisation tasks for the first call only
first_call: if (firstcall) then first_call: if (firstcall) then
@ -52,7 +54,6 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
! Set a simple limiter following the boundary of the data grid ! Set a simple limiter following the boundary of the data grid
simple_limiter: block simple_limiter: block
use const_and_precisions, only : cm use const_and_precisions, only : cm
use types, only : contour
! Avoid clipping out the launcher ! Avoid clipping out the launcher
real(wp_) :: R_launch, R_max real(wp_) :: R_launch, R_max
@ -60,37 +61,37 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
R_max = max(R_launch, R(mr)) R_max = max(R_launch, R(mr))
! Convert to a closed contour ! Convert to a closed contour
data%equilibrium%limiter = contour( & limiter = contour(Rmin=params%misc%rwall, Rmax=R_max, &
Rmin=params%misc%rwall, Rmax=R_max, zmin=z(1), zmax=z(mz)) zmin=z(1), zmax=z(mz))
end block simple_limiter end block simple_limiter
end if first_call end if first_call
! Set MHD equilibrium data ! Set MHD equilibrium data
init_equilibrium: block init_equilibrium: block
use equilibrium, only : set_equil_spline use gray_equil, only : eqdsk_data
use types, only : contour
! Copy argument arrays type(eqdsk_data) :: data
data%equilibrium%rv = r
data%equilibrium%zv = z ! Build EQDSK structure from in-memory data
data%equilibrium%rax = rax data%grid_r = r
data%equilibrium%rvac = rax data%grid_z = z
data%equilibrium%zax = zax data%axis = [rax, zax]
data%equilibrium%psinr = psrad data%r_ref = rax
data%equilibrium%fpol = fpol data%psi = psrad
data%equilibrium%psia = psia data%fpol = fpol
data%equilibrium%psin = psin data%psi_a = psia
data%equilibrium%qpsi = qpsi data%psi_map = psin
data%equilibrium%boundary = contour(rbnd, zbnd) data%q = qpsi
data%boundary = contour(rbnd, zbnd)
! Compute splines ! Compute splines
call set_equil_spline(params%equilibrium, data%equilibrium, err) call equil%init(params%equilibrium, data, err)
if (err /= 0) return if (err /= 0) return
end block init_equilibrium end block init_equilibrium
! Compute splines of kinetic profiles ! Compute splines of kinetic profiles
call plasma%init(params, psrad, te, dne, zeff, err) call plasma%init(params, equil, psrad, te, dne, zeff, err)
if (err /= 0) return if (err /= 0) return
! Set wave launcher parameters ! Set wave launcher parameters
@ -108,7 +109,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
end block init_antenna end block init_antenna
! Call main subroutine for the ibeam-th beam ! Call main subroutine for the ibeam-th beam
call gray_main(params, data, plasma, res, err, rhout=sqrt(psrad)) call gray_main(params, equil, plasma, limiter, res, err, rhout=sqrt(psrad))
write_debug_files: block write_debug_files: block
integer :: i, err integer :: i, err
@ -120,14 +121,6 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
end do end do
end block write_debug_files end block write_debug_files
! Free memory
free_memory: block
use equilibrium, only : unset_equil_spline
! Unset global variables of the `equilibrium` module
call unset_equil_spline
end block free_memory
! Copy over the results ! Copy over the results
pabs = res%pabs pabs = res%pabs
icd = res%icd icd = res%icd

30
src/gray_params.f90
View File

@ -167,30 +167,6 @@ module gray_params
integer, allocatable :: active_tables(:) ! IDs of output tables to fill in integer, allocatable :: active_tables(:) ! IDs of output tables to fill in
end type end type
! MHD equilibrium data
type equilibrium_data
real(wp_), allocatable :: rv(:) ! R of the uniform grid
real(wp_), allocatable :: zv(:) ! Z of the uniform grid
type(contour) :: limiter ! limiter contour (wall)
type(contour) :: boundary ! boundary contour (plasma)
real(wp_), allocatable :: fpol(:) ! Poloidal current function
real(wp_), allocatable :: qpsi(:) ! Safety factor on the flux grid
real(wp_), allocatable :: psin(:,:) ! Poloidal flux on a uniform grid
real(wp_), allocatable :: psinr(:) ! Poloidal flux
real(wp_) :: psia ! Poloidal flux at edge - flux at magnetic axis
real(wp_) :: rvac ! Reference R₀ (B = B₀R₀/R without the plasma)
real(wp_) :: rax ! R of the magnetic axis
real(wp_) :: zax ! Z of the magnetic axis
end type
! Kinetic plasma profiles data
type profiles_data
real(wp_), allocatable :: psrad(:) ! Radial coordinate
real(wp_), allocatable :: terad(:) ! Electron temperature profile
real(wp_), allocatable :: derad(:) ! Electron density profile
real(wp_), allocatable :: zfc(:) ! Effective charge profile
end type
! All GRAY parameters ! All GRAY parameters
type gray_parameters type gray_parameters
type(antenna_parameters) :: antenna type(antenna_parameters) :: antenna
@ -221,12 +197,6 @@ module gray_params
"misc.rwall" & "misc.rwall" &
] ]
! All GRAY input data
type gray_data
type(equilibrium_data) :: equilibrium
type(profiles_data) :: profiles
end type
! Wrapper type for array of pointers ! Wrapper type for array of pointers
type table_ptr type table_ptr
type(table), pointer :: ptr => null() type(table), pointer :: ptr => null()

18
src/gray_plasma.f90
View File

@ -262,17 +262,18 @@ contains
end function numeric_zeff end function numeric_zeff
subroutine load_plasma(params, plasma, err) subroutine load_plasma(params, equil, plasma, err)
! Loads a generic plasma description from file (params%filenm)
use gray_params, only : gray_parameters use gray_params, only : gray_parameters
use gray_params, only : PROF_ANALYTIC, PROF_NUMERIC use gray_params, only : PROF_ANALYTIC, PROF_NUMERIC
use gray_params, only : RHO_TOR, RHO_POL, RHO_PSI use gray_params, only : RHO_TOR, RHO_POL, RHO_PSI
use gray_params, only : SCALE_COLLISION, SCALE_GREENWALD, SCALE_OFF use gray_params, only : SCALE_COLLISION, SCALE_GREENWALD, SCALE_OFF
use gray_equil, only : abstract_equil
use logger, only : log_error, log_debug use logger, only : log_error, log_debug
use equilibrium, only : frhopol
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(inout) :: params type(gray_parameters), intent(inout) :: params
class(abstract_equil), intent(in) :: equil
class(abstract_plasma), allocatable, intent(out) :: plasma class(abstract_plasma), allocatable, intent(out) :: plasma
integer, intent(out) :: err integer, intent(out) :: err
@ -362,14 +363,14 @@ contains
! Convert psi to ψ_n (normalised poloidal flux) ! Convert psi to ψ_n (normalised poloidal flux)
select case (params%profiles%irho) select case (params%profiles%irho)
case (RHO_TOR) ! psi is ρ_t = √Φ_n (toroidal flux) case (RHO_TOR) ! psi is ρ_t = √Φ_n (toroidal flux)
psi = [(frhopol(psi(i))**2, i=1,nrows)] psi = [(equil%tor2pol(psi(i))**2, i=1,nrows)]
case (RHO_POL) ! psi is ρ_p = √ψ_n (poloidal flux) case (RHO_POL) ! psi is ρ_p = √ψ_n (poloidal flux)
psi = psi**2 psi = psi**2
case (RHO_PSI) ! psi is already ψ_n case (RHO_PSI) ! psi is already ψ_n
end select end select
! Interpolate ! Interpolate
call np%init(params, psi, temp, dens, zeff, err) call np%init(params, equil, psi, temp, dens, zeff, err)
end block end block
allocate(plasma, source=np) allocate(plasma, source=np)
@ -381,17 +382,19 @@ contains
end subroutine load_plasma end subroutine load_plasma
subroutine init_splines(self, params, psi, temp, dens, zeff, err) subroutine init_splines(self, params, equil, psi, temp, dens, zeff, err)
! Computes splines for the plasma profiles data and stores them ! Computes splines for the plasma profiles data and stores them
! in their respective global variables, see the top of this file. ! in their respective global variables, see the top of this file.
! !
! `err` is 1 if I/O errors occured, 2 if other initialisation failed. ! `err` is 1 if I/O errors occured, 2 if other initialisation failed.
use gray_params, only : gray_parameters use gray_params, only : gray_parameters
use gray_equil, only : abstract_equil
use logger, only : log_debug, log_info, log_warning, log_error use logger, only : log_debug, log_info, log_warning, log_error
! subroutine arguments ! subroutine arguments
class(numeric_plasma), intent(out) :: self class(numeric_plasma), intent(out) :: self
type(gray_parameters), intent(inout) :: params type(gray_parameters), intent(inout) :: params
class(abstract_equil), intent(in) :: equil
real(wp_), dimension(:), intent(in) :: psi, temp, dens, zeff real(wp_), dimension(:), intent(in) :: psi, temp, dens, zeff
integer, intent(out) :: err integer, intent(out) :: err
@ -481,7 +484,6 @@ contains
! Note: if it does, the initial wave conditions become ! Note: if it does, the initial wave conditions become
! invalid as they are given assuming a vacuum (N=1) ! invalid as they are given assuming a vacuum (N=1)
block block
use equilibrium, only : pol_flux
use const_and_precisions, only : cm use const_and_precisions, only : cm
real(wp_) :: R, Z, psi real(wp_) :: R, Z, psi
@ -493,7 +495,7 @@ contains
! Get the poloidal flux at the launcher ! Get the poloidal flux at the launcher
! Note: this returns -1 when the data is not available ! Note: this returns -1 when the data is not available
call pol_flux(R, z, psi) call equil%pol_flux(R, z, psi)
if (psi > self%tail%start .and. psi < self%tail%end) then if (psi > self%tail%start .and. psi < self%tail%end) then
! Fall back to the midpoint of ψ₀ and the launcher ψ ! Fall back to the midpoint of ψ₀ and the launcher ψ

81
src/gray_tables.f90
View File

@ -110,16 +110,17 @@ contains
end subroutine init_tables end subroutine init_tables
function kinetic_profiles(params, plasma) result(tbl) function kinetic_profiles(params, equil, plasma) result(tbl)
! Generates the plasma kinetic profiles ! Generates the plasma kinetic profiles
use gray_params, only : gray_parameters, EQ_VACUUM use gray_params, only : gray_parameters, EQ_VACUUM
use gray_equil, only : abstract_equil
use gray_plasma, only : abstract_plasma use gray_plasma, only : abstract_plasma
use equilibrium, only : fq, frhotor
use magsurf_data, only : npsi, vajphiav use magsurf_data, only : npsi, vajphiav
! function arguments ! function arguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
class(abstract_plasma), intent(in) :: plasma class(abstract_plasma), intent(in) :: plasma
type(table) :: tbl type(table) :: tbl
@ -145,7 +146,7 @@ contains
do i = 0, npsi + ntail do i = 0, npsi + ntail
rho_p = i * drho_p rho_p = i * drho_p
rho_t = frhotor(rho_p) rho_t = equil%pol2tor(rho_p)
psi_n = rho_p**2 psi_n = rho_p**2
if (psi_n < 1) then if (psi_n < 1) then
J_phi = vajphiav(i+1) * 1.e-6_wp_ J_phi = vajphiav(i+1) * 1.e-6_wp_
@ -153,24 +154,26 @@ contains
J_phi = 0 J_phi = 0
end if end if
call plasma%density(psi_n, dens, ddens) call plasma%density(psi_n, dens, ddens)
call tbl%append([psi_n, rho_t, dens, plasma%temp(psi_n), fq(psi_n), J_phi]) call tbl%append([psi_n, rho_t, dens, plasma%temp(psi_n), &
equil%safety(psi_n), J_phi])
end do end do
end function kinetic_profiles end function kinetic_profiles
function ec_resonance(params, B_res) result(tbl) function ec_resonance(params, equil, B_res) result(tbl)
! Generates the EC resonance surface table ! Generates the EC resonance surface table
use const_and_precisions, only : comp_eps use const_and_precisions, only : comp_eps
use gray_params, only : gray_parameters, EQ_VACUUM use gray_params, only : gray_parameters, EQ_VACUUM
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield use gray_equil, only : abstract_equil
use magsurf_data, only : npsi use magsurf_data, only : npsi
use types, only : item use types, only : item
use cniteq, only : cniteq_f use cniteq, only : cniteq_f
! function arguments ! function arguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
real(wp_), intent(in) :: B_res real(wp_), intent(in) :: B_res
type(table) :: tbl type(table) :: tbl
@ -190,11 +193,11 @@ contains
real(wp_), dimension(npsi) :: R, z real(wp_), dimension(npsi) :: R, z
! Build a regular (R, z) grid ! Build a regular (R, z) grid
dr = (rmxm - rmnm - comp_eps)/(npsi - 1) dr = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(npsi - 1)
dz = (zmxm - zmnm)/(npsi - 1) dz = (equil%z_range(2) - equil%z_range(1))/(npsi - 1)
do j=1,npsi do j=1,npsi
R(j) = comp_eps + rmnm + dr*(j - 1) R(j) = comp_eps + equil%r_range(1) + dr*(j - 1)
z(j) = zmnm + dz*(j - 1) z(j) = equil%z_range(1) + dz*(j - 1)
end do end do
! B_tot on psi grid ! B_tot on psi grid
@ -202,7 +205,7 @@ contains
B_min = +1.0e30_wp_ B_min = +1.0e30_wp_
do k = 1, npsi do k = 1, npsi
do j = 1, npsi do j = 1, npsi
call bfield(R(j), z(k), B_R, B_z, B_phi) call equil%b_field(R(j), z(k), B_R, B_z, B_phi)
B_tot(j,k) = sqrt(B_R**2 + B_z**2 + B_phi**2) B_tot(j,k) = sqrt(B_R**2 + B_z**2 + B_phi**2)
if(B_tot(j,k) >= B_max) B_max = B_tot(j,k) if(B_tot(j,k) >= B_max) B_max = B_tot(j,k)
if(B_tot(j,k) <= B_min) B_min = B_tot(j,k) if(B_tot(j,k) <= B_min) B_min = B_tot(j,k)
@ -227,17 +230,18 @@ contains
end function ec_resonance end function ec_resonance
function inputs_maps(params, plasma, B_res, X_norm) result(tbl) function inputs_maps(params, equil, plasma, B_res, X_norm) result(tbl)
! Generates 2D maps of several input quantities ! Generates 2D maps of several input quantities
use const_and_precisions, only : comp_eps, cm, degree use const_and_precisions, only : comp_eps, cm, degree
use gray_params, only : gray_parameters, EQ_VACUUM use gray_params, only : gray_parameters, EQ_VACUUM
use gray_equil, only : abstract_equil
use gray_plasma, only : abstract_plasma use gray_plasma, only : abstract_plasma
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield
use magsurf_data, only : npsi use magsurf_data, only : npsi
! function arguments ! function arguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
class(abstract_plasma), intent(in) :: plasma class(abstract_plasma), intent(in) :: plasma
real(wp_), intent(in) :: B_res ! resonant magnetic field, e/m_eω real(wp_), intent(in) :: B_res ! resonant magnetic field, e/m_eω
real(wp_), intent(in) :: X_norm ! X normalisation, e²/ε₀m_eω² real(wp_), intent(in) :: X_norm ! X normalisation, e²/ε₀m_eω²
@ -261,18 +265,18 @@ contains
Npl0 = sin(params%antenna%beta*degree) ! initial value of N∥ Npl0 = sin(params%antenna%beta*degree) ! initial value of N∥
! Build a regular (R, z) grid ! Build a regular (R, z) grid
dR = (rmxm - rmnm - comp_eps)/(npsi - 1) dR = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(npsi - 1)
dz = (zmxm - zmnm)/(npsi - 1) dz = (equil%z_range(2) - equil%z_range(1))/(npsi - 1)
do j = 1, npsi do j = 1, npsi
R(j) = comp_eps + rmnm + dR*(j - 1) R(j) = comp_eps + equil%r_range(1) + dR*(j - 1)
z(j) = zmnm + dz*(j - 1) z(j) = equil%z_range(1) + dz*(j - 1)
end do end do
do j = 1, npsi do j = 1, npsi
Npl = Npl0 * R0/r(j) Npl = Npl0 * R0/r(j)
do k = 1, npsi do k = 1, npsi
call pol_flux(r(j), z(k), psi_n) call equil%pol_flux(r(j), z(k), psi_n)
call bfield(r(j), z(k), B_R, B_z, B_phi) call equil%b_field(r(j), z(k), B_R, B_z, B_phi)
call plasma%density(psi_n, ne, dne) call plasma%density(psi_n, ne, dne)
B = sqrt(B_R**2 + B_phi**2 + B_z**2) B = sqrt(B_R**2 + B_phi**2 + B_z**2)
X = X_norm*ne X = X_norm*ne
@ -286,13 +290,13 @@ contains
end function inputs_maps end function inputs_maps
subroutine find_flux_surfaces(qvals, params, tbl) subroutine find_flux_surfaces(qvals, params, equil, tbl)
! Finds the ψ for a set of values of q and stores the ! Finds the ψ for a set of values of q and stores the
! associated surfaces to the flux surface table ! associated surfaces to the flux surface table
use gray_params, only : gray_parameters use gray_params, only : gray_parameters
use equilibrium, only : fq, frhotor, rmaxis, zmaxis, zbsup, zbinf use gray_equil, only : abstract_equil
use magsurf_data, only : contours_psi, npoints use magsurf_data, only : npoints
use logger, only : log_info use logger, only : log_info
use minpack, only : hybrj1 use minpack, only : hybrj1
use types, only : table use types, only : table
@ -300,6 +304,7 @@ contains
! subroutine arguments ! subroutine arguments
real(wp_), intent(in) :: qvals(:) real(wp_), intent(in) :: qvals(:)
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
type(table), intent(inout) :: tbl type(table), intent(inout) :: tbl
! local variables ! local variables
@ -319,7 +324,7 @@ contains
! searching near the boundary in case q is not monotonic. ! searching near the boundary in case q is not monotonic.
sol = [0.8_wp_] ! first guess sol = [0.8_wp_] ! first guess
! Solve fq(ψ_n) = qvals(i) for ψ_n ! Solve q(ψ_n) = qvals(i) for ψ_n
call hybrj1(equation, n=1, x=sol, fvec=fvec, fjac=fjac, & call hybrj1(equation, n=1, x=sol, fvec=fvec, fjac=fjac, &
ldfjac=1, tol=1e-3_wp_, info=info, wa=wa, lwa=7) ldfjac=1, tol=1e-3_wp_, info=info, wa=wa, lwa=7)
! Solution ! Solution
@ -334,19 +339,20 @@ contains
real(wp_), dimension(npoints) :: R_cont, z_cont real(wp_), dimension(npoints) :: R_cont, z_cont
real(wp_) :: R_hi, R_lo, z_hi, z_lo real(wp_) :: R_hi, R_lo, z_hi, z_lo
! Guesses for the high,low horzizontal-tangent points ! Guesses for the high,low horizontal-tangent points
R_hi = rmaxis; R_hi = equil%axis(1)
z_hi = (zbsup + zmaxis)/2 z_hi = (equil%z_boundary(2) + equil%axis(2))/2
R_lo = rmaxis R_lo = equil%axis(1)
z_lo = (zbinf + zmaxis)/2 z_lo = (equil%z_boundary(1) + equil%axis(2))/2
call contours_psi(params, psi_n, R_cont, z_cont, R_hi, z_hi, R_lo, z_lo) call equil%flux_contour(psi_n, params%misc%rwall, &
R_cont, z_cont, R_hi, z_hi, R_lo, z_lo)
call store_flux_surface(tbl, psi_n, R_cont, z_cont) call store_flux_surface(tbl, psi_n, R_cont, z_cont)
end block end block
! Log the values found for ψ_n, ρ_p, ρ_t ! Log the values found for ψ_n, ρ_p, ρ_t
rho_p = sqrt(psi_n) rho_p = sqrt(psi_n)
rho_t = frhotor(rho_p) rho_t = equil%pol2tor(rho_p)
write (msg, '(4(x,a,"=",g0.3))') & write (msg, '(4(x,a,"=",g0.3))') &
'q', qvals(i), 'ψ_n', psi_n, 'ρ_p', rho_p, 'ρ_t', rho_t 'q', qvals(i), 'ψ_n', psi_n, 'ρ_p', rho_p, 'ρ_t', rho_t
call log_info(msg, mod='gray_tables', proc='find_flux_surfaces') call log_info(msg, mod='gray_tables', proc='find_flux_surfaces')
@ -368,10 +374,10 @@ contains
if (flag == 1) then if (flag == 1) then
! return f(x) ! return f(x)
f(1) = fq(x(1)) - qvals(i) f(1) = equil%safety(x(1)) - qvals(i)
else else
! return f'(x), computed numerically ! return f'(x), computed numerically
df(1,1) = (fq(x(1) + e) - fq(x(1) - e)) / (2*e) df(1,1) = (equil%safety(x(1) + e) - equil%safety(x(1) - e)) / (2*e)
end if end if
end subroutine end subroutine
@ -402,7 +408,7 @@ contains
end subroutine store_flux_surface end subroutine store_flux_surface
subroutine store_ray_data(params, tables, & subroutine store_ray_data(params, equil, tables, &
i, jk, s, P0, pos, psi_n, B, b_n, k0, & i, jk, s, P0, pos, psi_n, B, b_n, k0, &
N_pl, N_pr, N, N_pr_im, n_e, T_e, & N_pl, N_pr, N, N_pr_im, n_e, T_e, &
alpha, tau, dIds, nhm, nhf, iokhawa, & alpha, tau, dIds, nhm, nhf, iokhawa, &
@ -410,15 +416,16 @@ contains
! Stores some ray variables and local quantities ! Stores some ray variables and local quantities
use const_and_precisions, only : degree, cm use const_and_precisions, only : degree, cm
use equilibrium, only : frhotor use gray_equil, only : abstract_equil
use gray_params, only : gray_parameters, gray_tables use gray_params, only : gray_parameters, gray_tables
use beamdata, only : unfold_index use beamdata, only : unfold_index
! subroutine arguments ! subroutine arguments
! tables where to store the data ! tables where to store the data
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
type(gray_tables), intent(inout) :: tables class(abstract_equil), intent(in) :: equil
type(gray_tables), intent(inout) :: tables
integer, intent(in) :: i, jk ! step, ray index integer, intent(in) :: i, jk ! step, ray index
real(wp_), intent(in) :: s ! arclength real(wp_), intent(in) :: s ! arclength
@ -453,7 +460,7 @@ contains
if (jk == 1 .and. tables%central_ray%active) then if (jk == 1 .and. tables%central_ray%active) then
phi_deg = atan2(pos_m(2), pos_m(1)) / degree phi_deg = atan2(pos_m(2), pos_m(1)) / degree
if(psi_n >= 0 .and. psi_n <= 1) then if(psi_n >= 0 .and. psi_n <= 1) then
rho_t = frhotor(sqrt(psi_n)) rho_t = equil%pol2tor(sqrt(psi_n))
else else
rho_t = 1.0_wp_ rho_t = 1.0_wp_
end if end if

157
src/magsurf_data.f90
View File

@ -95,16 +95,16 @@ contains
end subroutine dealloc_surfvec end subroutine dealloc_surfvec
subroutine compute_flux_averages(params, tables) subroutine compute_flux_averages(params, equil, tables)
use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv
use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge,zbsup,zbinf, & use dierckx, only : regrid, coeff_parder
bfield,frhotor,fq,tor_curr,psia,pol_flux
use dierckx, only : regrid,coeff_parder
use types, only : table, wrap use types, only : table, wrap
use gray_params, only : gray_parameters, gray_tables use gray_params, only : gray_parameters, gray_tables
use gray_equil, only : abstract_equil
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
type(gray_tables), intent(inout), optional :: tables type(gray_tables), intent(inout), optional :: tables
! local constants ! local constants
@ -130,8 +130,7 @@ contains
real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam
real(wp_), dimension(lwrk) :: wrk real(wp_), dimension(lwrk) :: wrk
real(wp_), dimension(:), allocatable :: rctemp,zctemp real(wp_), dimension(:), allocatable :: rctemp,zctemp
! common/external functions/variables real(wp_) :: fpolv, ddpsidrr, ddpsidzz
real(wp_) :: fpolv,ddpsidrr,ddpsidzz
npsi=nnintp npsi=nnintp
npoints = 2*ncnt+1 npoints = 2*ncnt+1
@ -162,30 +161,31 @@ contains
dffhlam(nlam)=-99999.0_wp_ dffhlam(nlam)=-99999.0_wp_
jp=1 jp=1
anorm=2.0_wp_*pi*rmaxis/abs(btaxis) anorm=2.0_wp_*pi*equil%axis(1)/abs(equil%b_axis)
dvdpsi=2.0_wp_*pi*anorm dvdpsi=2.0_wp_*pi*anorm
dadpsi=2.0_wp_*pi/abs(btaxis) dadpsi=2.0_wp_*pi/abs(equil%b_axis)
b2av=btaxis**2 b2av=equil%b_axis**2
ratio_cdator=abs(btaxis/btrcen) ratio_cdator=abs(equil%b_axis/equil%b_centre)
ratio_cdbtor=1.0_wp_ ratio_cdbtor=1.0_wp_
ratio_pltor=1.0_wp_ ratio_pltor=1.0_wp_
fc=1.0_wp_ fc=1.0_wp_
call pol_flux(rmaxis, zmaxis, ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz) call equil%pol_flux(equil%axis(1), equil%axis(2), &
qq=abs(btaxis)/sqrt(ddpsidrr*ddpsidzz)/psia ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz)
ajphiav=-ccj*(ddpsidrr+ddpsidzz)*psia/rmaxis qq=abs(equil%b_axis)/sqrt(ddpsidrr*ddpsidzz)/equil%psi_a
ajphiav=-ccj*(ddpsidrr+ddpsidzz)*equil%psi_a/equil%axis(1)
psicon(1)=0.0_wp_ psicon(1)=0.0_wp_
rcon(:,1)=rmaxis rcon(:,1)=equil%axis(1)
zcon(:,1)=zmaxis zcon(:,1)=equil%axis(2)
pstab(1)=0.0_wp_ pstab(1)=0.0_wp_
rpstab(1)=0.0_wp_ rpstab(1)=0.0_wp_
vcurrp(1)=0.0_wp_ vcurrp(1)=0.0_wp_
vajphiav(1)=ajphiav vajphiav(1)=ajphiav
bmxpsi(1)=abs(btaxis) bmxpsi(1)=abs(equil%b_axis)
bmnpsi(1)=abs(btaxis) bmnpsi(1)=abs(equil%b_axis)
bav(1)=abs(btaxis) bav(1)=abs(equil%b_axis)
rbav(1)=1.0_wp_ rbav(1)=1.0_wp_
rri(1)=rmaxis rri(1)=equil%axis(1)
varea(1)=0.0_wp_ varea(1)=0.0_wp_
vvol(1)=0.0_wp_ vvol(1)=0.0_wp_
vratjpl(1)=ratio_pltor vratjpl(1)=ratio_pltor
@ -196,21 +196,22 @@ contains
dadrhotv(1)=0.0_wp_ dadrhotv(1)=0.0_wp_
dvdrhotv(1)=0.0_wp_ dvdrhotv(1)=0.0_wp_
rup=rmaxis rup=equil%axis(1)
rlw=rmaxis rlw=equil%axis(1)
zup=zmaxis+(zbsup-zmaxis)/10.0_wp_ zup=equil%axis(2)+(equil%z_boundary(2)-equil%axis(2))/10.0_wp_
zlw=zmaxis-(zmaxis-zbinf)/10.0_wp_ zlw=equil%axis(2)-(equil%axis(2)-equil%z_boundary(1))/10.0_wp_
do jp=2,npsi do jp=2,npsi
height=dble(jp-1)/dble(npsi-1) height=dble(jp-1)/dble(npsi-1)
if(jp.eq.npsi) height=0.9999_wp_ if(jp.eq.npsi) height=0.9999_wp_
rhopjp=height rhopjp=height
psinjp=height*height psinjp=height*height
rhotjp=frhotor(rhopjp) rhotjp=equil%pol2tor(rhopjp)
if (rhotjp /= rhotjp) print *, 'Nan!!!!!!!!!!' if (rhotjp /= rhotjp) print *, 'Nan!!!!!!!!!!'
psicon(jp)=height psicon(jp)=height
call contours_psi(params, psinjp, rctemp, zctemp, rup, zup, rlw, zlw) call equil%flux_contour(psinjp, params%misc%rwall, &
rctemp, zctemp, rup, zup, rlw, zlw)
rcon(:,jp) = rctemp rcon(:,jp) = rctemp
zcon(:,jp) = zctemp zcon(:,jp) = zctemp
@ -226,8 +227,8 @@ contains
bmmx=-1.0e+30_wp_ bmmx=-1.0e+30_wp_
bmmn=1.0e+30_wp_ bmmn=1.0e+30_wp_
ajphi0 = tor_curr(rctemp(1),zctemp(1)) ajphi0 = equil%tor_curr(rctemp(1), zctemp(1))
call bfield(rctemp(1), zctemp(1), brr, bzz, bphi) call equil%b_field(rctemp(1), zctemp(1), brr, bzz, bphi)
fpolv=bphi*rctemp(1) fpolv=bphi*rctemp(1)
btot0=sqrt(bphi**2+brr**2+bzz**2) btot0=sqrt(bphi**2+brr**2+bzz**2)
bpoloid0=sqrt(brr**2+bzz**2) bpoloid0=sqrt(brr**2+bzz**2)
@ -237,8 +238,8 @@ contains
do inc=1,npoints-1 do inc=1,npoints-1
inc1=inc+1 inc1=inc+1
dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) dla=sqrt((rctemp(inc)-equil%axis(1))**2+(zctemp(inc)-equil%axis(2))**2)
dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) dlb=sqrt((rctemp(inc1)-equil%axis(1))**2+(zctemp(inc1)-equil%axis(2))**2)
dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+(zctemp(inc1)-zctemp(inc))**2) dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+(zctemp(inc1)-zctemp(inc))**2)
drc=(rctemp(inc1)-rctemp(inc)) drc=(rctemp(inc1)-rctemp(inc))
@ -253,8 +254,8 @@ contains
! compute line integrals on the contour psi=psinjp=height^2 ! compute line integrals on the contour psi=psinjp=height^2
rpsim=rctemp(inc1) rpsim=rctemp(inc1)
zpsim=zctemp(inc1) zpsim=zctemp(inc1)
call bfield(rpsim, zpsim, brr, bzz) call equil%b_field(rpsim, zpsim, brr, bzz)
ajphi = tor_curr(rpsim,zpsim) ajphi = equil%tor_curr(rpsim, zpsim)
bphi=fpolv/rpsim bphi=fpolv/rpsim
btot=sqrt(bphi**2+brr**2+bzz**2) btot=sqrt(bphi**2+brr**2+bzz**2)
bpoloid=sqrt(brr**2+bzz**2) bpoloid=sqrt(brr**2+bzz**2)
@ -310,7 +311,7 @@ contains
bmxpsi(jp)=bmmx bmxpsi(jp)=bmmx
bmnpsi(jp)=bmmn bmnpsi(jp)=bmmn
rri(jp)=bav(jp)/abs(fpolv*r2iav) rri(jp)=bav(jp)/abs(fpolv*r2iav)
ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen)) ratio_cdator=abs(b2av*riav/(fpolv*r2iav*equil%b_centre))
ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav)) ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav))
ratio_pltor=abs(bbav*riav/(fpolv*r2iav)) ratio_pltor=abs(bbav*riav/(fpolv*r2iav))
vratjpl(jp)=ratio_pltor vratjpl(jp)=ratio_pltor
@ -318,8 +319,8 @@ contains
vratjb(jp)=ratio_cdbtor vratjb(jp)=ratio_cdbtor
qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi)) qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi))
qqv(jp)=qq qqv(jp)=qq
dadrhotv(jp)=phitedge*rhotjp/fq(psinjp)*dadpsi/pi dadrhotv(jp)=equil%phi_a*rhotjp/equil%safety(psinjp)*dadpsi/pi
dvdrhotv(jp)=phitedge*rhotjp/fq(psinjp)*dvdpsi/pi dvdrhotv(jp)=equil%phi_a*rhotjp/equil%safety(psinjp)*dvdpsi/pi
! computation of fraction of circulating/trapped fraction fc, ft ! computation of fraction of circulating/trapped fraction fc, ft
! and of function H(lambda,rhop) ! and of function H(lambda,rhop)
@ -397,8 +398,8 @@ contains
do jp=1,npsi do jp=1,npsi
if (.not. tables%flux_averages%active) exit if (.not. tables%flux_averages%active) exit
call tables%flux_averages%append([ & call tables%flux_averages%append([ &
rpstab(jp), frhotor(rpstab(jp)), bav(jp), bmxpsi(jp), & rpstab(jp), equil%pol2tor(rpstab(jp)), bav(jp), bmxpsi(jp), &
bmnpsi(jp), varea(jp), vvol(jp), vcurrp(jp), vajphiav(jp), & bmnpsi(jp), varea(jp), vvol(jp), vcurrp(jp), vajphiav(jp), &
ffc(jp), vratja(jp), vratjb(jp), qqv(jp)]) ffc(jp), vratja(jp), vratjb(jp), qqv(jp)])
end do end do
@ -444,88 +445,4 @@ contains
end subroutine fluxval end subroutine fluxval
subroutine contours_psi(params, h, rcn, zcn, rup, zup, rlw, zlw)
use const_and_precisions, only : wp_, pi
use gray_params, only : gray_parameters
use logger, only : log_warning
use dierckx, only : profil, sproota
use equilibrium, only : model, frhotor, psi_spline, &
kspl, find_htg_point
! local constants
integer, parameter :: mest=4
! subroutine arguments
type(gray_parameters), intent(in) :: params
real(wp_), intent(in) :: h
real(wp_), intent(out) :: rcn(:), zcn(:)
real(wp_), intent(inout) :: rup, zup, rlw, zlw
! local variables
integer :: npoints,np,ic,ier,iopt,m
real(wp_) :: ra,rb,za,zb,th,zc
real(wp_), dimension(mest) :: zeroc
real(wp_), dimension(psi_spline%nknots_x) :: czc
npoints=size(rcn)
np=(npoints-1)/2
th=pi/dble(np)
if (params%equilibrium%iequil<2) then
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
rb=rlw
za=zup
zb=zlw
call find_htg_point(ra,za,rup,zup,h)
call find_htg_point(rb,zb,rlw,zlw,h)
rcn(1)=rlw
zcn(1)=zlw
rcn(npoints)=rlw
zcn(npoints)=zlw
rcn(np+1)=rup
zcn(np+1)=zup
do ic=2,np
zc=zlw+(zup-zlw)*(1.0_wp_-cos(th*(ic-1)))/2.0_wp_
iopt=1
call profil(iopt, psi_spline%knots_x, psi_spline%nknots_x, &
psi_spline%knots_y, psi_spline%nknots_y, &
psi_spline%coeffs, kspl, kspl, zc, &
psi_spline%nknots_x, czc, ier)
if (ier > 0) then
block
character(256) :: msg
write(msg, '(a, a, g0)') &
'when computing ψ(R,z) contour `profil` returned ier=', ier
call log_warning(msg, mod='magsurf_data', proc='contours_psi')
end block
end if
call sproota(h, psi_spline%knots_x, psi_spline%nknots_x, &
czc, zeroc, mest, m, ier)
if (zeroc(1) > params%misc%rwall) then
rcn(ic)=zeroc(1)
rcn(npoints+1-ic)=zeroc(2)
else
rcn(ic)=zeroc(2)
rcn(npoints+1-ic)=zeroc(3)
end if
zcn(ic)=zc
zcn(npoints+1-ic)=zc
end do
end if
end subroutine contours_psi
end module magsurf_data end module magsurf_data

188
src/main.f90
View File

@ -2,12 +2,14 @@ program main
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
use logger, only : INFO, ERROR, WARNING, set_log_level, log_message use logger, only : INFO, ERROR, WARNING, set_log_level, log_message
use utils, only : dirname use utils, only : dirname
use types, only : contour
use gray_cli, only : cli_options, parse_cli_options, & use gray_cli, only : cli_options, parse_cli_options, &
parse_param_overrides parse_param_overrides
use gray_core, only : gray_main use gray_core, only : gray_main
use gray_equil, only : abstract_equil, load_equil
use gray_plasma, only : abstract_plasma, load_plasma use gray_plasma, only : abstract_plasma, load_plasma
use gray_params, only : gray_parameters, gray_data, gray_results, & use gray_params, only : gray_parameters, gray_results, &
read_gray_params, read_gray_config, & read_gray_params, read_gray_config, &
make_gray_header make_gray_header
implicit none implicit none
@ -18,8 +20,9 @@ program main
! gray_main subroutine arguments ! gray_main subroutine arguments
type(gray_parameters) :: params ! Inputs type(gray_parameters) :: params ! Inputs
type(gray_data) :: data ! class(abstract_equil), allocatable :: equil !
class(abstract_plasma), allocatable :: plasma ! class(abstract_plasma), allocatable :: plasma !
type(contour) :: limiter !
type(gray_results) :: results ! Outputs type(gray_results) :: results ! Outputs
integer :: err ! Exit code integer :: err ! Exit code
@ -68,7 +71,7 @@ program main
associate (p => params%raytracing) associate (p => params%raytracing)
if (p%nrayr < 5) then if (p%nrayr < 5) then
p%igrad = 0 p%igrad = 0
call log_message(level=WARNING, mod='main', proc='main', & call log_message(level=WARNING, mod='main', &
msg='nrayr < 5 ⇒ optical case only') msg='nrayr < 5 ⇒ optical case only')
end if end if
if (p%nrayr == 1) p%nrayth = 1 if (p%nrayr == 1) p%nrayth = 1
@ -78,29 +81,33 @@ program main
end if end if
end associate end associate
! Read the input data and set the global variables ! Read and initialise the equilibrium and limiter objects
! of the respective module. Note: order matters. call load_equil(params%equilibrium, equil, limiter, err)
call init_equilibrium(params, data, err)
if (err /= 0) call exit(err) if (err /= 0) call exit(err)
! Read and initialise the plasma state object ! Read and initialise the plasma state object
call load_plasma(params, plasma, err) call load_plasma(params, equil, plasma, err)
if (err /= 0) call exit(err) if (err /= 0) call exit(err)
call init_misc(params, data) ! Create a simple limiter, if necessary
if (.not. allocated(limiter%R) .or. params%raytracing%ipass < 0) then
call log_message('using fallback limiter', level=INFO, mod='main')
params%raytracing%ipass = abs(params%raytracing%ipass)
limiter = make_fallback_limiter(params, equil)
end if
! Get back to the initial directory ! Get back to the initial directory
err = chdir(cwd) err = chdir(cwd)
if (allocated(opts%sum_filelist)) then if (allocated(opts%sum_filelist)) then
call log_message(level=INFO, mod='main', msg='summing profiles') call log_message(level=INFO, mod='main', msg='summing profiles')
call sum_profiles(params, opts%sum_filelist, opts%output_dir, results) call sum_profiles(params, equil, opts%sum_filelist, opts%output_dir, results)
else else
! Activate the given output tables ! Activate the given output tables
params%misc%active_tables = opts%tables params%misc%active_tables = opts%tables
! Finally run the simulation ! Finally run the simulation
call gray_main(params, data, plasma, results, err) call gray_main(params, equil, plasma, limiter, results, err)
end if end if
print_res: block print_res: block
@ -136,68 +143,10 @@ program main
end do end do
end block write_res end block write_res
! Free memory
free_mem: block
use equilibrium, only : unset_equil_spline
call unset_equil_spline
end block free_mem
end block no_save end block no_save
contains contains
subroutine init_equilibrium(params, data, err)
! Reads the MHD equilibrium file (either in the G-EQDSK format
! or an analytical description) and initialises the respective
! GRAY parameters and data.
use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
use equilibrium, only : read_equil_an, read_eqdsk, change_cocos, &
set_equil_an, set_equil_spline, scale_equil
use logger, only : log_debug
! subroutine arguments
type(gray_parameters), intent(inout) :: params
type(gray_data), intent(out) :: data
integer, intent(out) :: err
select case (params%equilibrium%iequil)
case (EQ_VACUUM)
call log_debug('vacuum, no MHD equilibrium', &
mod='main', proc='init_equilibrium')
case (EQ_ANALYTICAL)
call log_debug('loading analytical file', &
mod='main', proc='init_equilibrium')
call read_equil_an(params%equilibrium%filenm, &
params%raytracing%ipass, &
data%equilibrium, err)
if (err /= 0) return
case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
call log_debug('loading G-EQDK file', &
mod='main', proc='init_equilibrium')
call read_eqdsk(params%equilibrium, data%equilibrium, err)
if (err /= 0) return
call change_cocos(data%equilibrium, params%equilibrium%icocos, 3)
end select
! Rescale B, I and/or force their signs
call scale_equil(params%equilibrium, data%equilibrium)
! Set global variables
select case (params%equilibrium%iequil)
case (EQ_ANALYTICAL)
call set_equil_an
case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
call log_debug('computing splines...', mod='main', proc='init_equilibrium')
call set_equil_spline(params%equilibrium, data%equilibrium, err)
if (err /= 0) return
call log_debug('splines computed', mod='main', proc='init_equilibrium')
end select
end subroutine init_equilibrium
subroutine init_antenna(params, err) subroutine init_antenna(params, err)
! Reads the wave launcher file (containing the wave frequency, launcher ! Reads the wave launcher file (containing the wave frequency, launcher
! position, direction and beam description) and initialises the respective ! position, direction and beam description) and initialises the respective
@ -227,82 +176,61 @@ contains
end subroutine init_antenna end subroutine init_antenna
subroutine init_misc(params, data) function make_fallback_limiter(params, equil) result(limiter)
! Performs miscellanous initial tasks, before the gray_main subroutine. ! Creates a simple rectangular limiter:
use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, & !
EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL ! (R_wall, z_launch+Δz)╔═════╗(R_max, z_launch+Δz)
use types, only : contour ! ║4 3║
use const_and_precisions, only : cm, comp_huge ! ║ ║
! ║1 2║
! (R_wall, z_launch-Δz)╚═════╝(R_max, z_launch-Δz)
!
! Note: R_max = { +∞ for vacuum,
! { R₀+a for analytical model,
! { grid_r[-1] for numerical data.
!
use logger, only : log_info
use const_and_precisions, only : cm
! subroutine arguments ! function arguments
type(gray_parameters), intent(inout) :: params type(gray_parameters), intent(in) :: params
type(gray_data), intent(inout) :: data class(abstract_equil), intent(in) :: equil
type(contour) :: limiter
! Build a basic limiter if one is not provided by equilibrium.txt ! local variables
if (.not. allocated(data%equilibrium%limiter%R) & real(wp_) :: R_launch, z_launch, R_max, delta_z
.or. params%raytracing%ipass < 0) then
! Restore sign of ipass ! Launcher coordinates
params%raytracing%ipass = abs(params%raytracing%ipass) R_launch = norm2(params%antenna%pos(1:2)) * cm
z_launch = params%antenna%pos(3) * cm
block ! Max vertical distance a ray can travel
real(wp_) :: R_launch, z_launch, R_max, delta_z delta_z = abs(params%raytracing%ipass) * &
params%raytracing%dst * params%raytracing%nstep * cm
! Launcher coordinates ! Avoid clipping out the launcher
R_launch = norm2(params%antenna%pos(1:2)) * cm R_max = max(R_launch, equil%r_range(2))
z_launch = params%antenna%pos(3) * cm
! Max vertical distance a ray can travel ! Build a rectangle
delta_z = abs(params%raytracing%ipass) * & limiter = contour( &
params%raytracing%dst * params%raytracing%nstep * cm Rmin=params%misc%rwall, Rmax=R_max, &
zmin=z_launch - delta_z, zmax=z_launch + delta_z)
! Max radius, either due to the plasma extent or equilibrium grid end function make_fallback_limiter
select case (params%equilibrium%iequil)
case (EQ_VACUUM)
! Use a very large R, ~ unbounded
R_max = comp_huge
case (EQ_ANALYTICAL)
! use R₀+a
block
use equilibrium, only : model
R_max = model%R0 + model%a
end block
case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
! use max R of the grid
R_max = data%equilibrium%rv(size(data%equilibrium%rv))
end select
! Avoid clipping out the launcher
R_max = max(R_launch, R_max)
! Convert to a list of R,z:
!
! (R_wall, z_launch+Δz)╔═════╗(R_max, z_launch+Δz)
! ║4 3║
! ║ ║
! ║1 2║
! (R_wall, z_launch-Δz)╚═════╝(R_max, z_launch-Δz)
!
data%equilibrium%limiter = contour( &
Rmin=params%misc%rwall, Rmax=R_max, &
zmin=z_launch - delta_z, zmax=z_launch + delta_z)
end block
end if
end subroutine init_misc
subroutine sum_profiles(params, filelist, output_dir, results) subroutine sum_profiles(params, equil, filelist, output_dir, results)
! Combines the EC profiles obtained from multiple gray runs ! Combines the EC profiles obtained from multiple gray runs
! (of different beams) and recomputes the summary table ! (of different beams) and recomputes the summary table
use gray_params, only : gray_parameters use gray_params, only : gray_parameters
use gray_equil, only : abstract_equil
use magsurf_data, only : compute_flux_averages, dealloc_surfvec use magsurf_data, only : compute_flux_averages, dealloc_surfvec
use pec, only : pec_init, postproc_profiles, dealloc_pec, & use pec, only : pec_init, postproc_profiles, dealloc_pec, &
rhot_tab rhot_tab
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(inout) :: params type(gray_parameters), intent(inout) :: params
class(abstract_equil), intent(in) :: equil
character(len=*), intent(in) :: filelist, output_dir character(len=*), intent(in) :: filelist, output_dir
type(gray_results), intent(inout) :: results type(gray_results), intent(inout) :: results
@ -322,10 +250,10 @@ contains
integer, allocatable :: beams(:) integer, allocatable :: beams(:)
! Initialise the magsurf_data module ! Initialise the magsurf_data module
call compute_flux_averages(params) call compute_flux_averages(params, equil)
! Initialise the output profiles ! Initialise the output profiles
call pec_init(params%output) call pec_init(params%output, equil)
associate(nrho =>params%output%nrho) associate(nrho =>params%output%nrho)
allocate(jphi(nrho), currins(nrho), pins(nrho), rtin(nrho), rpin(nrho)) allocate(jphi(nrho), currins(nrho), pins(nrho), rtin(nrho), rpin(nrho))
@ -419,7 +347,7 @@ contains
results%icd = currins(params%output%nrho) results%icd = currins(params%output%nrho)
! Recompute the summary values ! Recompute the summary values
call postproc_profiles( & call postproc_profiles(equil, &
results%Pabs, results%Icd, rhot_tab, results%dPdV, jphi, & results%Pabs, results%Icd, rhot_tab, results%dPdV, jphi, &
rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip, & rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip, &
rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx) rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx)

78
src/multipass.f90
View File

@ -2,26 +2,28 @@ module multipass
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
use polarization, only : pol_limit, field_to_ellipse use polarization, only : pol_limit, field_to_ellipse
use reflections, only : wall_refl use reflections, only : wall_refl
use equilibrium, only : bfield use gray_equil, only : abstract_equil
implicit none implicit none
contains contains
subroutine plasma_in(i, x, N, Bres, sox, cpl, psi, chi, iop, ext, eyt, perfect) subroutine plasma_in(i, equil, x, N, Bres, sox, cpl, &
psi, chi, iop, ext, eyt, perfect)
! Computes the ray polarisation and power couplings when it enteres the plasma ! Computes the ray polarisation and power couplings when it enteres the plasma
use const_and_precisions, only : cm use const_and_precisions, only : cm
! subroutine arguments ! subroutine arguments
integer, intent(in) :: i ! ray index integer, intent(in) :: i ! ray index
real(wp_), intent(in) :: x(3), N(3) ! position, refactive index class(abstract_equil), intent(in) :: equil ! MHD equilibrium
real(wp_), intent(in) :: Bres ! resonant B field real(wp_), intent(in) :: x(3), N(3) ! position, refactive index
integer, intent(in) :: sox ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X real(wp_), intent(in) :: Bres ! resonant B field
real(wp_), intent(out) :: cpl(2) ! power coupling vector (O, X) integer, intent(in) :: sox ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
real(wp_), intent(out) :: psi, chi ! polarisation ellipse angles real(wp_), intent(out) :: cpl(2) ! power coupling vector (O, X)
integer, intent(inout) :: iop(:) ! inside/outside plasma flag real(wp_), intent(out) :: psi, chi ! polarisation ellipse angles
complex(wp_), intent(inout) :: ext(:), eyt(:) ! ray polarisation vector (e_x, e_y) integer, intent(inout) :: iop(:) ! inside/outside plasma flag
logical, intent(in) :: perfect ! whether to assume perfect coupling complex(wp_), intent(inout) :: ext(:), eyt(:) ! ray polarisation vector (e_x, e_y)
logical, intent(in) :: perfect ! whether to assume perfect coupling
! local variables ! local variables
real(wp_) :: R, z, cosphi, sinphi, B_phi, B_R, B_z real(wp_) :: R, z, cosphi, sinphi, B_phi, B_R, B_z
@ -37,7 +39,7 @@ contains
z = x(3) * cm z = x(3) * cm
cosphi = x(1)/R * cm cosphi = x(1)/R * cm
sinphi = x(2)/R * cm sinphi = x(2)/R * cm
call bfield(R, z, B_R, B_z, B_phi) call equil%b_field(R, z, B_R, B_z, B_phi)
B(1) = B_R*cosphi - B_phi*sinphi B(1) = B_R*cosphi - B_phi*sinphi
B(2) = B_R*sinphi + B_phi*cosphi B(2) = B_R*sinphi + B_phi*cosphi
B(3) = B_z B(3) = B_z
@ -70,16 +72,17 @@ contains
end subroutine plasma_in end subroutine plasma_in
subroutine plasma_out(i, xv, anv, bres, sox, iop, ext, eyt) subroutine plasma_out(i, equil, xv, anv, bres, sox, iop, ext, eyt)
! Ray exits plasma ! Ray exits plasma
! subroutine arguments ! subroutine arguments
integer, intent(in) :: i ! ray index integer, intent(in) :: i ! ray index
real(wp_), intent(in) :: xv(3), anv(3) class(abstract_equil), intent(in) :: equil ! MHD equilibrium
real(wp_), intent(in) :: bres real(wp_), intent(in) :: xv(3), anv(3)
integer, intent(in) :: sox real(wp_), intent(in) :: bres
integer, intent(inout) :: iop(:) ! in/out plasma flag integer, intent(in) :: sox
complex(wp_), intent(out) :: ext(:), eyt(:) integer, intent(inout) :: iop(:) ! in/out plasma flag
complex(wp_), intent(out) :: ext(:), eyt(:)
! local variables ! local variables
real(wp_) :: rm, csphi, snphi, bphi, br, bz real(wp_) :: rm, csphi, snphi, bphi, br, bz
@ -91,7 +94,7 @@ contains
rm=sqrt(xmv(1)**2+xmv(2)**2) rm=sqrt(xmv(1)**2+xmv(2)**2)
csphi=xmv(1)/rm csphi=xmv(1)/rm
snphi=xmv(2)/rm snphi=xmv(2)/rm
call bfield(rm,xmv(3),br,bz,bphi) call equil%b_field(rm,xmv(3),br,bz,bphi)
bv(1)=br*csphi-bphi*snphi bv(1)=br*csphi-bphi*snphi
bv(2)=br*snphi+bphi*csphi bv(2)=br*snphi+bphi*csphi
bv(3)=bz bv(3)=bz
@ -99,34 +102,35 @@ contains
end subroutine plasma_out end subroutine plasma_out
subroutine wall_out(i, wall, ins, xv, anv, dst, bres, sox, & subroutine wall_out(i, equil, wall, ins, xv, anv, dst, bres, &
psipol1, chipol1, iow, iop, ext, eyt) sox, psipol1, chipol1, iow, iop, ext, eyt)
! Ray exits vessel ! Ray exits vessel
use types, only : contour use types, only : contour
! subroutine arguments ! subroutine arguments
integer, intent(in) :: i ! ray index integer, intent(in) :: i ! ray index
type(contour), intent(in) :: wall ! wall contour class(abstract_equil), intent(in) :: equil ! MHD equilibrium
logical, intent(in) :: ins ! inside plasma? (ins=1 plasma/wall overlap) type(contour), intent(in) :: wall ! wall contour
real(wp_), intent(inout) :: xv(3), anv(3) logical, intent(in) :: ins ! inside plasma? (ins=1 plasma/wall overlap)
real(wp_), intent(in) :: dst ! step size real(wp_), intent(inout) :: xv(3), anv(3)
real(wp_), intent(in) :: bres real(wp_), intent(in) :: dst ! step size
integer, intent(in) :: sox real(wp_), intent(in) :: bres
real(wp_), intent(out) :: psipol1, chipol1 integer, intent(in) :: sox
integer, intent(inout) :: iow(:) ! in/out vessel and plasma flags real(wp_), intent(out) :: psipol1, chipol1
integer, intent(inout) :: iop(:) ! in/out vessel and plasma flags integer, intent(inout) :: iow(:) ! in/out vessel and plasma flags
complex(wp_), intent(inout) :: ext(:), eyt(:) integer, intent(inout) :: iop(:) ! in/out vessel and plasma flags
complex(wp_), intent(inout) :: ext(:), eyt(:)
! local variables ! local variables
integer :: irfl integer :: irfl
real(wp_), dimension(3) :: xvrfl,anvrfl,walln real(wp_), dimension(3) :: xvrfl,anvrfl,walln
complex(wp_) :: ext1,eyt1 complex(wp_) :: ext1,eyt1
iow(i)=iow(i)+1 ! out->in iow(i) = iow(i) + 1 ! out->in
if(ins) call plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! plasma-wall overlapping if(ins) call plasma_out(i, equil, xv, anv, bres, sox, iop, ext, eyt) ! plasma-wall overlapping
call wall_refl(wall, xv-dst*anv, anv, ext(i), eyt(i), & call wall_refl(wall, xv-dst*anv, anv, ext(i), eyt(i), &
xvrfl, anvrfl, ext1, eyt1, walln, irfl) ! ray reflects at wall xvrfl, anvrfl, ext1, eyt1, walln, irfl) ! ray reflects at wall
ext(i) = ext1 ! save parameters at wall reflection ext(i) = ext1 ! save parameters at wall reflection
eyt(i) = eyt1 eyt(i) = eyt1
xv = xvrfl xv = xvrfl
anv = anvrfl anv = anvrfl

38
src/pec.f90
View File

@ -10,13 +10,14 @@ module pec
contains contains
subroutine pec_init(params, rt_in) subroutine pec_init(params, equil, rt_in)
use gray_params, only : output_parameters, RHO_POL use gray_params, only : output_parameters, RHO_POL
use equilibrium, only : frhotor, frhopol use gray_equil, only : abstract_equil
use magsurf_data, only : fluxval use magsurf_data, only : fluxval
! subroutine arguments ! subroutine arguments
type(output_parameters), intent(in) :: params type(output_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
real(wp_), intent(in), optional :: rt_in(params%nrho) real(wp_), intent(in), optional :: rt_in(params%nrho)
! local variables ! local variables
@ -60,12 +61,12 @@ contains
end if end if
if (params%ipec == RHO_POL) then if (params%ipec == RHO_POL) then
rhop_tab(it) = rt rhop_tab(it) = rt
rhot_tab(it) = frhotor(rt) rhot_tab(it) = equil%pol2tor(rt)
rhop1 = rt1 rhop1 = rt1
else else
rhot_tab(it) = rt rhot_tab(it) = rt
rhop_tab(it) = frhopol(rt) rhop_tab(it) = equil%tor2pol(rt)
rhop1 = frhopol(rt1) rhop1 = equil%tor2pol(rt1)
end if end if
! psi grid at mid points, size n+1, for use in pec_tab ! psi grid at mid points, size n+1, for use in pec_tab
@ -242,23 +243,24 @@ contains
end subroutine pec_tab end subroutine pec_tab
subroutine postproc_profiles(pabs,currt,rhot_tab,dpdv,ajphiv, & subroutine postproc_profiles(equil, pabs, currt, rhot_tab, dpdv, ajphiv, &
rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip, & rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip, &
rhotp, drhotp, rhotjfi, drhotjfi, dpdvmx,ajmxfi, ratjamx, ratjbmx) rhotp, drhotp, rhotjfi, drhotjfi, dpdvmx,ajmxfi, ratjamx, ratjbmx)
! Radial average values over power and current density profile ! Radial average values over power and current density profile
use const_and_precisions, only : pi, comp_eps use const_and_precisions, only : pi, comp_eps
use equilibrium, only : frhopol use gray_equil, only : abstract_equil
use magsurf_data, only : fluxval use magsurf_data, only : fluxval
real(wp_), intent(in) :: pabs,currt class(abstract_equil), intent(in) :: equil
real(wp_), intent(in) :: rhot_tab(:) real(wp_), intent(in) :: pabs,currt
real(wp_), intent(in) :: dpdv(:), ajphiv(:) real(wp_), intent(in) :: rhot_tab(:)
real(wp_), intent(out) :: rhotpav, drhotpav, dpdvp real(wp_), intent(in) :: dpdv(:), ajphiv(:)
real(wp_), intent(out) :: rhotjava, drhotjava, ajphip real(wp_), intent(out) :: rhotpav, drhotpav, dpdvp
real(wp_), intent(out) :: rhotp, drhotp, dpdvmx real(wp_), intent(out) :: rhotjava, drhotjava, ajphip
real(wp_), intent(out) :: rhotjfi, drhotjfi, ajmxfi real(wp_), intent(out) :: rhotp, drhotp, dpdvmx
real(wp_), intent(out) :: ratjamx, ratjbmx real(wp_), intent(out) :: rhotjfi, drhotjfi, ajmxfi
real(wp_), intent(out) :: ratjamx, ratjbmx
real(wp_) :: sccsa, ratjplmx, rhopjava, rhoppav real(wp_) :: sccsa, ratjplmx, rhopjava, rhoppav
real(wp_) :: rhotjav, rhot2pav, rhot2java, dvdrhotav, dadrhotava real(wp_) :: rhotjav, rhot2pav, rhot2java, dvdrhotav, dadrhotava
@ -287,8 +289,8 @@ contains
drhotpav = sqrt(8._wp_*(max(rhot2pav -rhotpav**2,comp_eps))) drhotpav = sqrt(8._wp_*(max(rhot2pav -rhotpav**2,comp_eps)))
drhotjava = sqrt(8._wp_*(max(rhot2java-rhotjava**2,comp_eps))) drhotjava = sqrt(8._wp_*(max(rhot2java-rhotjava**2,comp_eps)))
rhoppav = frhopol(rhotpav) rhoppav = equil%tor2pol(rhotpav)
rhopjava = frhopol(rhotjava) rhopjava = equil%tor2pol(rhotjava)
if (pabs > 0) then if (pabs > 0) then
call fluxval(rhoppav,dvdrhot=dvdrhotav) call fluxval(rhoppav,dvdrhot=dvdrhotav)

18
src/splines.f90
View File

@ -58,6 +58,7 @@ module splines
procedure :: init_deriv => spline_2d_init_deriv procedure :: init_deriv => spline_2d_init_deriv
procedure :: deriv => spline_2d_deriv procedure :: deriv => spline_2d_deriv
procedure :: transform => spline_2d_transform procedure :: transform => spline_2d_transform
final :: spline_2d_final
end type end type
! A simple 1D linear interpolation ! A simple 1D linear interpolation
@ -507,6 +508,23 @@ contains
end subroutine spline_2d_deinit end subroutine spline_2d_deinit
subroutine spline_2d_final(self)
! Deallocates pointer components in a spline_2d
type(spline_2d), intent(inout) :: self
if (allocated(self%partial)) then
deallocate(self%partial(1, 0)%ptr)
deallocate(self%partial(0, 1)%ptr)
deallocate(self%partial(1, 1)%ptr)
deallocate(self%partial(2, 0)%ptr)
deallocate(self%partial(0, 2)%ptr)
deallocate(self%partial)
end if
end subroutine spline_2d_final
function spline_2d_eval(self, x, y) result(z) function spline_2d_eval(self, x, y) result(z)
! Evaluates the spline at (x, y) ! Evaluates the spline at (x, y)
use dierckx, only : fpbisp use dierckx, only : fpbisp