!> Deals with most types of equilibria to generate GS2's expected geometrical parameters !> Input: !> rhoc :: minor radius of surface of interest (see irho below) !> use_large_aspect :: default false. Choose true to get "s-alpha" equilibria (see local_eq) !> Equilibrium options : Choose 1 of: !> local_eq : Analytic parameterization of local equilibrium !> ppl_eq = .true. :: use PPL style NetCDF equilibrium (Menard) !> transp_eq = .true. :: use TRXPL style NetCDF equilibrium (Menard) !> gen_eq = .true. :: use GA style NetCDF equilibrium (TOQ) !> chs_eq = .true. :: Use CHEASE generated equilibrium !> efit_eq = .true. :: use EFIT equilibrium !> dfit_eq = .true. :: use dipole equilibrium !> -------------------> [note: if any of the above except local_eq are true set eqfile = input file name] !> irho :: choose flux surface label !> case (1) :: sqrt(toroidal flux)/sqrt(toroidal flux of LCFS) !> case (2) :: diameter/(diameter of LCFS). recommended !> case (3) :: poloidal flux/(poloidal flux of LCFS) !> case (4) :: rho_mid? Only used for the vacuum ring dipole, dfit_eq = T !> !> equal_arc :: .true. makes grad_par coefficient a constant. !> !> New feature: Given an actual equilibrium, allow s_hat and d pressure/d rho !> to vary self-consistently (a la Greene and Chance, Bishop). We use the !> Bishop formalism [C. M. Bishop, et al., NF, Vol. 24, 1579 (1984).] !> !> bishop :: uses Bishop relations to generate coefficients. !> case(0) :: do not use Bishop relations -- primarily for testing !> case(1) :: use Bishop relations with actual equilibrium shat, p' !> case(2) :: use Bishop relations with new shat, and p' from p_prime_input !> case(3) :: use Bishop relations with new shat, L_p from below !> s_hat_input :: new magnetic shear, rho/q*dq/drho !> invLp_input :: new -1/pressure * d pressure/d rho !> case(4) :: use Bishop relations with new shat, beta' from below !> s_hat_input :: new magnetic shear, rho/q*dq/drho !> beta_prime_input :: new d beta/d rho !> case(5) :: use Bishop relations with new shat, alpha from below !> s_hat_input :: new magnetic shear, rho/q*dq/drho !> alpha_input :: new alpha = -q**2 * rmaj * d pressure/d rho !> case(6) :: use Bishop relations with equilibrium shat, beta' from below !> beta_prime_input :: new d beta/d rho !> case(7) :: use Bishop relations with equilibrium shat, beta' from below !> dp_mult :: d pressure/d rho = equilibrium gradient * dp_mult !> (In each case, the definition of rho is determined by irho above) !> !> nperiod :: number of cells of width 2 pi in theta. !> !> geoType :: determines the analytic geometry specification !> case (0) :: the tilted Miller geometry: generalization of Miller et al. PoP 1998 !> (see Chapter 3.2.1 of J. Ball MIT Masters thesis) !> case (1) :: the global MHD solution geometry: calculated by inverting the large !> aspect ratio, constant-current MHD solution !> (see J. Ball "Optimized up-down asymmetry to drive fast intrinsic rotation in tokamak reactors") !> case (2) :: the generalized ellipticity geometry: the polar formula for an ellipse generalized !> to arbitrary shaping effect mode number (see Chapter 5.1.2 of J. Ball Oxford PhD thesis) !> case (3) :: the Fourier series geometry: Fourier decomposition of the minor radius as a function !> of poloidal angle (see Chapter 5.1.1 of J. Ball Oxford PhD thesis) !> rmaj :: major radius of magnetic axis (in units of aminor) if not local_eq !> or major radius of center of flux surface of interest (units of aminor) otherwise !> r_geo :: major radius of last closed flux surface or !> if local_eq: sets the reference magnetic field to be the value of the toroidal !> field at R=r_geo on the flux surface of interest !> aSurf :: the real-space flux surface label of the surface with deltam and deltan shaping as well as !> a major radius of Rmaj (only used when geoType={1}) !> shift :: d rmaj/drho. (Typically < 0.) (only used when geoType={0,2,3}) !> -> Alternative definitions: !> Raxis :: major radial location of the magnetic axis. (only used when geoType={1}) !> shiftVert :: d zmaj/drho. (< 0 for a vertical shift) (only used when geoType={0,2,3}) !> -> Alternative definitions: !> Zaxis :: axial location of the magnetic axis. (only used when geoType={1}) !> mMode :: the first poloidal mode number (only used when geoType={2,3}) !> nMode :: the second poloidal mode number (only used when geoType={2,3}) !> deltam :: magnitude of the first poloidal shaping effect (only used when geoType={2,3}) !> -> Alternative definitions: !> delta2 :: identical to akappa (only used when geoType={1}) !> akappa :: flux surface elongation (only used when geoType={0}) !> deltan :: magnitude of the second poloidal shaping effect (only used when geoType={2,3}) !> -> Alternative definitions: !> delta3 :: magnitude of the triangulaity-like shaping effect (only used when geoType={1}) !> tri :: flux surface triangulaity (only used when geoType={0}) !> deltampri :: d(deltam)/d(rhoc) (only used when geoType={2,3}) !> -> Alternative definitions: !> akappapri :: d(akappa)/d(rhoc) (only used when geoType={0}) !> deltanpri :: d(deltan)/d(rhoc) (only used when geoType={2,3}) !> -> Alternative definitions: !> tripri :: d(tri)/d(rhoc) (only used when geoType={0}) !> thetam :: tilt angle of the first poloidal shaping effect (only used when geoType={2,3}) !> -> Alternative definitions: !> theta2 :: tilt angle of the flux surface elongation (only used when geoType={1}) !> thetak :: tilt angle of the flux surface elongation (only used when geoType={0}) !> thetan :: tilt angle of the first poloidal shaping effect (only used when geoType={2,3}) !> -> Alternative definitions: !> theta3 :: tilt angle of the triangularity-like shaping effect (only used when geoType={1}) !> thetad :: tilt angle of the flux surface triangularity (only used when geoType={0}) !> qinp :: safety factor q at surface of interest (local_eq = T only) !> shat :: d q/drho at surface of interest !> !> delrho :: numerical parameter. Should be "small enough". Typically 0.001 ok. !> !> force_sym :: true -> assume up-down symmetric equilibrium. !> writelots :: .true. for more output !> !> Use: !> Include this module, set input variables appropriately and call eikcoefs. !> !> I do not recommend calling other routines, but a few are available for !> diagnostic purposes. module geometry use geo_utils, only: abstract_geo_type, EQFILE_LENGTH, flux_surface_type implicit none private public :: eikcoefs, pbarfun, pbarofrho, rpofrho, diameter, run_eikcoefs_and_resample public :: eqfile, eqnormfile, irho, surf, verb, big public :: ppl_eq, transp_eq, chs_eq, efit_eq, gs2d_eq, idfit_eq, dfit_eq, local_eq, gen_eq public :: force_sym, writelots, equal_arc, bishop, use_large_aspect public :: p_prime_input, invLp_input, beta_prime_input, alpha_input, dp_mult public :: geom, eikcoefs_output_type, finish_geometry class(abstract_geo_type), allocatable :: geom ! Start of inputs type(flux_surface_type) :: surf real :: p_prime_input, invLp_input, beta_prime_input, alpha_input, dp_mult integer :: irho, bishop logical :: gen_eq, efit_eq, ppl_eq, local_eq, chs_eq, dfit_eq, idfit_eq, gs2d_eq, transp_eq logical :: writelots, equal_arc, force_sym, use_large_aspect character(len = EQFILE_LENGTH) :: eqfile, eqnormfile ! End of inputs real :: rpmin, rpmax integer :: big = 8 !> The single-period theta domain "half-grid" size. That is, the number of points in !> \((0, \pi]\). `theta(-nth:nth)` is the single central period \([-\pi, \pi]\). For !> EFIT and local equilbria, this is set from !> [[theta_grid_parameters_config_type:ntheta]]; for all other numerical equilbria, it !> is determined by the data file. Always even. integer :: nth !> The full theta domain "half-grid" size. That is, the total size of the !> [[geometry:theta]] grid is `(2 * ntgrid) + 1`, with array bounds !> `theta(-ntgrid:ntgrid)`. Equal to `(2*nperiod - 1)*nth` integer :: ntgrid !> Verbosity of print statements | Should we use runtime_tests:verbosity? integer :: verb = 2 !> If debug then print out lots of debug info. !> debug is true if (verb >= 3) logical :: debug = .false. !> A type to hold a collection of geometrical quantities as calculated by !> [[eikcoefs]] type eikcoefs_output_type integer :: ntheta real :: drhodpsi real :: kxfac real :: qsf real :: surfarea real :: dvdrhon real :: shat real :: dbetadrho real :: rhoc !> theta :: the theta grid of the results. The final grid has !> pi - 2*pi*nperiod <= theta <= 2*pi*nperiod - pi real, dimension(:), allocatable :: theta !> bmag :: |B| real, dimension(:), allocatable :: bmag !> gradpar :: coefficient of b_hat.grad operator real, dimension(:), allocatable :: gradpar !> grho :: grad(rho) real, dimension(:), allocatable :: grho !> cdrift :: part of coriolis drift operator independent of theta_0 real, dimension(:), allocatable :: cdrift !> cvdrift :: coefficient of v_par**2/Omega b_hat x (b_hat.grad b_hat).grad real, dimension(:), allocatable :: cvdrift !> gbdrift :: coefficient of mu/Omega b_hat x (grad B).grad operator real, dimension(:), allocatable :: gbdrift !> cdrift0 :: part of coriolis drift operator proportional to theta_0 real, dimension(:), allocatable :: cdrift0 !> cvdrift0 :: part of total curvature drift operator proportional to theta_0 real, dimension(:), allocatable :: cvdrift0 !> gbdrift0 :: part of total grad B drift operator proportional to theta_0 real, dimension(:), allocatable :: gbdrift0 !> gds2 :: part of grad_perp**2 operator independent of theta_0 real, dimension(:), allocatable :: gds2 !> gds21 :: part of grad_perp**2 operator proportional to theta_0 real, dimension(:), allocatable :: gds21 !> gds22 :: part of grad_perp**2 operator proportional to theta_0**2 real, dimension(:), allocatable :: gds22 !> gds23 :: part of v_E . grad theta operator independent of theta_0 real, dimension(:), allocatable :: gds23 !> gds24 :: part of v_E . grad theta operator proportional to theta_0 real, dimension(:), allocatable :: gds24 !> gds24_noq :: gds24/dqdrp real, dimension(:), allocatable :: gds24_noq !> rplot :: R(theta) - major radius on theta grid real, dimension(:), allocatable :: rplot !> zplot :: Z(theta) - height on theta grid real, dimension(:), allocatable :: zplot !> aplot :: a(theta) - alpha-phi on theta grid real, dimension(:), allocatable :: aplot !> rprime :: derivative wrt rho of R(theta) - major radius on theta grid real, dimension(:), allocatable :: rprime !> zprime :: derivative wrt rho Z(theta) - height on theta grid real, dimension(:), allocatable :: zprime !> aprime :: derivative wrt rho a(theta) - alpha-phi on theta grid real, dimension(:), allocatable :: aprime !> bpol :: Poloidal magnetic field strength real, dimension(:), allocatable :: bpol end type eikcoefs_output_type contains !> Run eikcoefs with `ntheta = ntheta_input` and then resample onto !> a theta grid with resolution `ntheta_target`. subroutine run_eikcoefs_and_resample(ntheta_input, ntheta_target, nperiod, results, job_id) implicit none !> This is the number of theta grid points per two pi period used in the !> geometry calculation (actually one less than this). integer, intent(in out) :: ntheta_input !> This is the desired number of theta grid points per two pi period !> in the resulting geometrical coefficients used elsewhere. integer, intent(in) :: ntheta_target integer, intent(in) :: nperiod type(eikcoefs_output_type), intent(out) :: results integer, intent(in), optional :: job_id type(eikcoefs_output_type) :: full_results real, dimension(:), allocatable :: input_index_space, target_index_space integer :: ig, ntgrid_target, ntgrid_input real :: spacing ! First run eikcoefs using requested resolution call eikcoefs(ntheta_input, nperiod, full_results, job_id) ntgrid_input = (size(full_results%theta) - 1) / 2 ntgrid_target =(ntheta_target*(2*nperiod-1)) / 2 ! Set scalar outputs results%ntheta = ntheta_target results%drhodpsi = full_results%drhodpsi results%kxfac = full_results%kxfac results%qsf = full_results%qsf results%surfarea = full_results%surfarea results%dvdrhon = full_results%dvdrhon results%shat = full_results%shat results%dbetadrho = full_results%dbetadrho results%rhoc = full_results%rhoc ! Make an index space array for the original theta ! grid. Spans [0,1] in size(theta)=2*ntgrid_input points allocate(input_index_space(-ntgrid_input:ntgrid_input)) spacing = 1.0 / (2 * ntgrid_input) do ig = -ntgrid_input, ntgrid_input input_index_space(ig) = (ig + ntgrid_input) * spacing end do ! Make new theta grid allocate(target_index_space(-ntgrid_target:ntgrid_target)) ! First make target index space grid spacing = 1.0 / (2 * ntgrid_target) do ig = -ntgrid_target, ntgrid_target target_index_space(ig) = (ig + ntgrid_target) * spacing end do ! Now interpolate theta grid onto target index space grid call resample( input_index_space, full_results%theta, target_index_space, results%theta) ! Now interpolate other output quantities onto target theta grid call resample(full_results%theta, full_results%bmag, results%theta, results%bmag) call resample(full_results%theta, full_results%gradpar, results%theta, results%gradpar) call resample(full_results%theta, full_results%grho, results%theta, results%grho) call resample(full_results%theta, full_results%cdrift, results%theta, results%cdrift) call resample(full_results%theta, full_results%cvdrift, results%theta, results%cvdrift) call resample(full_results%theta, full_results%gbdrift, results%theta, results%gbdrift) call resample(full_results%theta, full_results%cdrift0, results%theta, results%cdrift0) call resample(full_results%theta, full_results%cvdrift0, results%theta, results%cvdrift0) call resample(full_results%theta, full_results%gbdrift0, results%theta, results%gbdrift0) call resample(full_results%theta, full_results%gds2, results%theta, results%gds2) call resample(full_results%theta, full_results%gds21, results%theta, results%gds21) call resample(full_results%theta, full_results%gds22, results%theta, results%gds22) call resample(full_results%theta, full_results%gds23, results%theta, results%gds23) call resample(full_results%theta, full_results%gds24, results%theta, results%gds24) call resample(full_results%theta, full_results%gds24_noq, results%theta, results%gds24_noq) call resample(full_results%theta, full_results%rplot, results%theta, results%rplot) call resample(full_results%theta, full_results%zplot, results%theta, results%zplot) call resample(full_results%theta, full_results%aplot, results%theta, results%aplot) call resample(full_results%theta, full_results%rprime, results%theta, results%rprime) call resample(full_results%theta, full_results%zprime, results%theta, results%zprime) call resample(full_results%theta, full_results%aprime, results%theta, results%aprime) call resample(full_results%theta, full_results%bpol, results%theta, results%bpol) contains subroutine resample(x_orig, y_orig, x_new, y_new) use splines, only: new_spline, delete_spline, splint, spline real, dimension(:), intent(in) :: x_orig, y_orig real, dimension(-ntgrid_target:), intent(in) :: x_new real, dimension(:), allocatable, intent(out) :: y_new type(spline) :: the_spline integer :: i, lower, upper ! Spline the old data the_spline = new_spline(x_orig, y_orig) lower = lbound(x_new, 1) ; upper = ubound(x_new, 1) allocate(y_new(lower:upper)) do i = lower, upper y_new(i) = splint(x_new(i), the_spline) end do ! Get rid of the spline call delete_spline(the_spline) end subroutine resample end subroutine run_eikcoefs_and_resample !> FIXME : Add documentation subroutine eikcoefs (ntheta_input, nperiod, outputs, job_id_in) use constants, only: twopi use file_utils, only: open_output_file, close_output_file, finish_file_utils use file_utils, only: futils_initialized => initialized, init_file_utils use optionals, only: get_option_with_default use warning_helpers, only: is_zero implicit none !> Sets the number of theta grid points. If odd will be updated to be ntheta_input - 1 integer, intent(in out) :: ntheta_input integer, intent(in) :: nperiod !> A type containing all of the outputs of eikcoefs type(eikcoefs_output_type), intent(out) :: outputs !> job_id is the job id in a Trinity or list mode run. Only for debug integer, intent(in), optional :: job_id_in logical :: start_file_utils, dummy_flag integer :: job_id, i, iBmax, fort11_unit real, allocatable, dimension(:) :: seik, rgrid, invR, trip, dsdthet, rmajor, dummy real, allocatable, dimension(:,:) :: thgrad, rpgrad, crpgrad, grads real, allocatable, dimension(:,:) :: gradrptot, gradstot, gradztot, gradthtot, gradbtot real :: shat, dbetadrho, dqdrp, qval, dum, bi real :: rp, pbar, drhodrp, drhodrhod, drhodpsi, dpsidrp, dpsidrho ! Used to translate definition of the poloidal angle such that the location ! of the maximium magnetic field is at pi (see section 3.2.3 of Ball MIT ! Masters thesis or section 2.4 of "GS2 analytic geometry specification") real :: thetaShift job_id = get_option_with_default(job_id_in, -1) start_file_utils = .not. futils_initialized if (start_file_utils) then ! Initialize file utils with a run name we provide ! and by setting trin_run prevent it from searching ! the command line arguments. write(*,*) 'WARNING: Automatically initializing file_utils& & in eikcoefs. This is a hack to & & support older programs like eiktest. You should& & initialize file utils before calling eikcoefs' call init_file_utils(dummy_flag, name='geometry_module', trin_run=.true., input=.false.) end if !Create module level output file and get unit call open_output_file(fort11_unit,".fort.11") !Compute the initial constants debug = (verb > 2) if (debug) writelots = .true. if (debug) write(6,*) "eikcoefs: local_eq=",local_eq, 'jid=',job_id if (local_eq) then if (irho /= 2) write(*,*) 'Forcing irho = 2 for local_eq' irho = 2 if (bishop == 0) then write(*,*) 'Bishop = 0 not compatible with local_eq. Forcing bishop = 1' bishop = 1 !Not clear that this is a good choice for Miller end if end if if (debug) write(6,*) "eikcoefs: call check ",'jid=',job_id call check(gen_eq, efit_eq, ppl_eq, local_eq, dfit_eq, idfit_eq, chs_eq, gs2d_eq, transp_eq, bishop) call init_uniform_theta_grid(ntheta_input, outputs%theta, nperiod, ntgrid, nth) surf%nt = ntheta_input if (debug) write(6,*) 'eikcoefs: alloc_local_arrays' if (.not. allocated(seik)) call alloc_local_arrays(ntgrid) thetaShift = 0.0 allocate(outputs%bmag(-ntgrid:ntgrid), outputs%bpol(-ntgrid:ntgrid)) if (.not. allocated(geom)) call allocate_geom(geom) call initialise_geometry_and_get_bmag(geom, surf%r, thetaShift, outputs%theta, outputs%bmag, outputs%bpol) ! Start theta shift calc. This is used to modify the definition of ! theta such that the maxiumum magnetic field is at pi. This is ! necessary to correctly calculate bounce points. Can only shift with local_eq if (local_eq) then ! Locate the index of the maximum magnetic field in the range [-pi,pi) iBmax = maxloc(outputs%bmag(-nth:nth-1), dim = 1) - (nth + 1) ! If max is not at -pi and theta then shift theta. if (iBmax /= -nth) then thetaShift = outputs%theta(iBmax) - outputs%theta(-nth) call geom%finalise !Clean up geometry to allow reinitialisation write(*,*) "Bmax is not located at +/-pi, translating theta by",thetaShift," so that it is..." call initialise_geometry_and_get_bmag(geom, surf%r, thetaShift, outputs%theta, outputs%bmag, outputs%bpol) end if end if ! compute |grad rho|: allocate(outputs%grho(-ntgrid:ntgrid)) outputs%grho(-nth:nth) = hypot(rpgrad(-nth:nth,1), rpgrad(-nth:nth,2)) * abs(drhodrp) call periodic_copy(outputs%grho, 0., nperiod) if (force_sym) call sym(outputs%grho, 0, ntgrid) if (writelots) write(fort11_unit,*) 'dpsidrp= ',dpsidrp drhodpsi = drhodrp / dpsidrp ; dpsidrho = 1.0 / drhodpsi ! compute derivatives of eikonal, grad B and grad Z if (bishop /= 0) then call bishop_grads(surf%r, surf%rmaj, outputs%bmag, outputs%bpol) else call bishop_large_aspect_ratio_approximation(surf%r, surf%dr, ntheta_input) end if if (writelots) then write(fort11_unit,*) 'i,theta,grads1,2= ' do i = -ntgrid, ntgrid write(fort11_unit,*) i, outputs%theta(i), grads(i, 1), grads(i, 2) end do end if ! compute total grad S including toroidal part call gradstottgrid(invR, grads, gradstot) gradrptot(-nth:nth, 1:2) = rpgrad(-nth:nth, :) call periodic_copy(gradrptot(:, 1), 0., nperiod) call periodic_copy(gradrptot(:, 2), 0., nperiod) gradrptot(:, 3) = 0. ! gradthtot is grad theta since ballooning, theta goes from ! -ntheta*nperiod -> ntheta*nperiod, but grad theta is periodic in theta gradthtot(-nth:nth, 1:2) = thgrad(-nth:nth, :) call periodic_copy(gradthtot(:, 1), 0., nperiod) call periodic_copy(gradthtot(:, 2), 0., nperiod) gradthtot(:, 3) = 0. ! compute the magnitude squared of grad S allocate(outputs%gds2(-ntgrid:ntgrid), outputs%gds21(-ntgrid:ntgrid), & outputs%gds22(-ntgrid:ntgrid), outputs%gds23(-ntgrid:ntgrid), & outputs%gds24(-ntgrid:ntgrid), outputs%gds24_noq(-ntgrid:ntgrid)) call calculate_metric_terms(gradstot, gradrptot, gradthtot, dpsidrho, dqdrp, & qval, surf%r, outputs%bmag, force_sym, outputs%gds2, outputs%gds21, outputs%gds22, & outputs%gds23, outputs%gds24, outputs%gds24_noq) ! compute curvature, grad b, and coriolis drift terms allocate(outputs%gbdrift(-ntgrid:ntgrid), outputs%gbdrift0(-ntgrid:ntgrid)) allocate(outputs%cvdrift(-ntgrid:ntgrid), outputs%cvdrift0(-ntgrid:ntgrid)) allocate(outputs%cdrift(-ntgrid:ntgrid), outputs%cdrift0(-ntgrid:ntgrid)) call drift(rgrid, rp, gradstot, gradrptot, gradztot, gradbtot, dqdrp, dpsidrho, dpsidrp,& outputs%bmag, invR, outputs%theta, dbetadrho, bi, force_sym, nperiod, & outputs%gbdrift, outputs%gbdrift0, outputs%cvdrift, outputs%cvdrift0, & outputs%cdrift, outputs%cdrift0) !Recompute the magnetic field along theta if use_large_aspect - should this come earlier? if (use_large_aspect) then outputs%bmag(-nth:nth)= 1 / (1 + rgrid(-nth:nth) * cos(outputs%theta(-nth:nth)) / surf%rmaj) call periodic_copy(outputs%bmag, 0., nperiod) end if if (force_sym) call sym(outputs%bmag, 0, ntgrid) ! compute the grad-parallel coefficient allocate(outputs%gradpar(-ntgrid:ntgrid)) if (dfit_eq) then outputs%gradpar = trip / outputs%bmag else outputs%gradpar(-nth:nth) = -bi / (rmajor(-nth:nth)**2 * outputs%bmag(-nth:nth) * dsdthet(-nth:nth)) call periodic_copy(outputs%gradpar, 0., nperiod) end if if (force_sym) call sym(outputs%gradpar, 0, ntgrid) ! Note : We should call this before we potentially redefine the theta grid ! due to equal_arc below as some specific equilibrium types rely on theta ! having a consistent definition with that used to calculate rgrid initially. ! In particular this is important for Efit and GS2D equilibria (in eeq). Other ! equilibrium types use interpolation so shouldn't depend on the specific theta grid ! (although it's possible there could still be some interaction, given that one of the ! outputs of the below 'aplot' is given by `seik` which is evaluated on the original ! grid rather than the "new" one defined below). ! @note Check if restriction on cartesian geo types still holds. allocate(outputs%rplot(-ntgrid:ntgrid), outputs%rprime(-ntgrid:ntgrid)) allocate(outputs%zplot(-ntgrid:ntgrid), outputs%zprime(-ntgrid:ntgrid)) allocate(outputs%aplot(-ntgrid:ntgrid), outputs%aprime(-ntgrid:ntgrid)) call plotdata (rgrid, seik, grads, dpsidrho, surf%dr, outputs%theta, & outputs%rplot, outputs%zplot, outputs%aplot, & outputs%rprime, outputs%zprime, outputs%aprime) ! In the lines below we redefine theta to be theta', and gradpar ! to be gradpar'. This means that the values of theta in the gs2 ! output file are not values of the poloidal angle but values of theta_prime if (equal_arc) then allocate(dummy(-ntgrid:ntgrid)) call arclength (ntheta_input, nperiod, outputs%theta, outputs%gradpar, dummy) outputs%theta = dummy end if call calculate_surface_area(outputs%theta, outputs%gradpar, outputs%bmag, outputs%grho, & dpsidrho, surf%r, nth, outputs%surfarea, outputs%dvdrhon) if (writelots) write (fort11_unit,*) 'surfarea=', outputs%surfarea !Close output file call close_output_file(fort11_unit) outputs%kxfac = abs(qval * dpsidrho / surf%r) if (is_zero(qval) .and. (dfit_eq .or. idfit_eq)) then shat = max(1.e-7, shat) !RN> I think kxfac should include the sign. However, from GS2 convention, !RN> eventually kxfac = abs(kxfac). See get_grids in theta_grid. Even if the sign is !RN> effective, probably it has no effect due to symmetry. I just leave the sign here. outputs%kxfac = - dpsidrho outputs%gds22 = (outputs%grho * dpsidrho)**2 end if if (start_file_utils) call finish_file_utils ! Pack data into output type outputs%ntheta = ntheta_input ; outputs%drhodpsi = drhodpsi ; outputs%rhoc = surf%r outputs%qsf = qval ; outputs%shat = shat ; outputs%dbetadrho = dbetadrho contains !> Initialises the specific geometry backend subroutine initialise_geometry_and_get_bmag(geom, rhoc, thetaShift, theta, bmag, bpol) use geo_utils, only: geo_input_type class(abstract_geo_type), intent(in out) :: geom real, intent(in) :: rhoc, thetaShift real, dimension(-ntgrid:), intent(in) :: theta real, dimension(-ntgrid:), intent(out) :: bmag, bpol type(geo_input_type) :: geom_inputs real :: avgrmid, rbar, B_T0 integer :: fort78_unit, fort99_unit ! Initialise the specific geometry types (mostly reading data ! from file and calculating derived quantities). geom_inputs%surf = surf ; geom_inputs%thShift = thetaShift ; geom_inputs%big = big geom_inputs%eqfile = eqfile ; geom_inputs%eqnormfile = eqnormfile ; geom_inputs%theta = theta call geom%initialise(geom_inputs, rpmin, rpmax, surf%rmaj, B_T0, avgrmid) if (debug) write(6,*) "eikcoefs: rpmin,rpmax=", rpmin, rpmax if (.not. local_eq) surf%r_geo = geom%rcenter(rpmax) !Only used for writelots/checks if (surf%r_geo > surf%rmaj + 1.e-5) then write(6,*) write(6,*) 'Warning: the center of the LCFS is outboard of the ' write(6,*) 'center of the reference surface. This is probably wrong.' write(6,*) end if if (writelots) then write(fort11_unit,*) 'All lengths normalized to the avg midplane minor radius' if (irho /= 2) then write(fort11_unit,*) 'except rho, which always goes from 0 to 1.' write(fort11_unit,*) 'The rho normalization is determined by irho.' end if write(fort11_unit,*) 'avgrmid = ',avgrmid,' meters' if (.not. local_eq) write(fort11_unit,*) 'R_mag = ',surf%rmaj if (.not. local_eq) write(fort11_unit,*) 'B_T0 = ',B_T0 end if if (debug) write(6,*) "eikcoefs: find rp ", 'jid=',job_id rp = rpofrho(rhoc) ; drhodrp = drho_drp(rp, surf%dr) if (debug) write(6,*) "eikcoefs: rp=",rp, 'jid=',job_id if (debug) write(6,*) "eikcoefs: find rgrid" call rtgrid(rgrid, rp, theta) call periodic_copy(rgrid, 0.0, nperiod) pbar = pbarfun(rgrid(0), theta(0)) bi = geom%btori(pbar) ! Get the toroidal field function (i.e. fp = R B_T) if (writelots) then write(fort11_unit,*) 'rp of rho = ',rp write(fort11_unit,*) 'Rcenter(rho) = ',geom%rcenter(rp) write(fort11_unit,*) 'R_geo = ',surf%r_geo,' !or ',bi call open_output_file(fort78_unit,".fort.78") call open_output_file(fort99_unit,".fort.99") write(fort78_unit,*) '# theta_pol, R, Z' do i = -nth, nth write(fort78_unit,*) theta(i), geom%Rpos(rgrid(i), theta(i)), & geom%Zpos(rgrid(i), theta(i)) write(fort99_unit,*) 'i=', i, ' thet= ', theta(i), ' rgrid= ', rgrid(i) end do call close_output_file(fort78_unit) ; call close_output_file(fort99_unit) rbar = 0. do i = -nth, nth - 1 rbar = rbar + rgrid(i) * (theta(i + 1) - theta(i)) end do rbar = rbar / twopi write (*,*) 'r_bar = ', rbar write (*,*) write(fort11_unit,*) 'drhodrp=', drhodrp if (irho == 1) then drhodrhod = drho_drhod(rp, drhodrp, surf%dr) write(fort11_unit,*) 'drho_drhod=', drhodrhod end if end if if (debug) write(fort11_unit,*) 'eikcoefs: various logs', gen_eq, ppl_eq, chs_eq, efit_eq, dfit_eq, idfit_eq, gs2d_eq, transp_eq if (debug) write(6,*) "eikcoefs: call rmajortgrid" call rmajortgrid(rgrid, theta, rmajor) ; call periodic_copy(rmajor, 0., nperiod) if (use_large_aspect) then invR = 1.0 / surf%rmaj else call invRtgrid(rgrid, theta, invR) ; call periodic_copy(invR, 0., nperiod) end if ! compute gradient of rp in R, Z, phi (cartesian - crpgrad) ! Note that rpgrad is given in (R,Z,phi) coordinates if bishop == 0 ! and (rho,l,phi) coordinates otherwise ! Also, rpgrad is grad psi if not local_eq and grad rho if local_eq -- MAB ! When bishop > 0 ! rpgrad(:,1) is |grad rp| ; rpgrad(:,2) is 0 ! crpgrad(:,1) is drp/dR ; crpgrad(:,2) is drp/dZ ! ! bishop = 0 is high aspect ratio coefficients for debugging -- EGH call geom%gradient(rgrid, theta, crpgrad, 'P', dum, nth, ntgrid, .false.) call periodic_copy(crpgrad(:,1), 0., nperiod) call periodic_copy(crpgrad(:,2), 0., nperiod) ! compute coordinate gradients call geom%gradient(rgrid, theta, rpgrad, 'P', dum, nth, ntgrid, bishop /= 0) call geom%gradient(rgrid, theta, thgrad, 'T', dum, nth, ntgrid, bishop /= 0) call tripprod2dtgrid(rpgrad, thgrad, invR, trip) call periodic_copy(trip, 0., nperiod) ! compute eikonal S if (local_eq) qval = geom%qfun(pbar) !leq:qfun ignores input pbar if (debug) write(6,*) "eikcoefs: call eikonal" call eikonal(theta, invR, trip, qval, bi, seik, dsdthet, dpsidrp) if (debug) write(6,*) "eikcoefs: done eikonal" if (writelots) write(fort11_unit,*) 'q= ', qval if (local_eq) trip = trip * dpsidrp if (bishop > 0) then !Following bits _could_ be calculated for all bishop call B_pol(invR, crpgrad, dpsidrp, bpol, nth) !<Should this use rpgrad for local? call periodic_copy(bpol, 0., nperiod) call B_mod(invR, bpol, bmag, bi, nth) call periodic_copy(bmag, 0., nperiod) end if !Note bmag not set for bishop == 0 end subroutine initialise_geometry_and_get_bmag subroutine bishop_grads(rhoc, rmaj, bmag, bpol) real, intent(in) :: rhoc, rmaj real, dimension(-ntgrid:), intent(in) :: bmag, bpol logical, parameter :: save_bishop = .true. integer :: bishop_unit real :: a_b, b_b, c_b, dp, dp_new, di, di_new, s_hat, s_hat_new, pressure, tmp real, allocatable, dimension(:) :: a_bish, b_bish, c_bish, da_bish, db_bish, dc_bish, & zoftheta, dzdl, dsdl, ltheta, th_bish, invRpol allocate(a_bish(-ntgrid:ntgrid), b_bish(-ntgrid:ntgrid), c_bish(-ntgrid:ntgrid), & da_bish(-ntgrid:ntgrid), db_bish(-ntgrid:ntgrid), dc_bish(-ntgrid:ntgrid), & zoftheta(-ntgrid:ntgrid), dzdl(-ntgrid:ntgrid), dsdl(-ntgrid:ntgrid), & ltheta(-ntgrid:ntgrid), th_bish(-ntgrid:ntgrid), invRpol(-ntgrid:ntgrid)) call th_bishop(crpgrad, th_bish, nth) call periodic_copy(th_bish, 0., nperiod) call loftheta(rgrid, outputs%theta, ltheta) call periodic_copy(ltheta, ltheta(nth) - ltheta(-nth), nperiod) call inv_Rpol(outputs%theta, th_bish, ltheta, invRpol, nth) call periodic_copy(invRpol, 0., nperiod) if (debug) write(6,*) "eikcoefs: gradients" if (debug) write(6,*) 'drhodpsi', drhodpsi if (debug) write(6,*) 'bi', bi if (verb > 5) write(6,*) 'rmajor', rmajor(:) if (verb > 5) write(6,*) 'bpol', bpol(:) if (verb > 5) write(6,*) rpgrad(-nth:nth,1) if (verb > 5) write(6,*) rpgrad(-nth:nth,2) if (verb > 5) write(6,*) thgrad(-nth:nth,1) if (verb > 5) write(6,*) thgrad(-nth:nth,2) if (verb > 5) write(6,*) 1.0 / trip(-nth:nth) if (verb > 5) write(6,*) -1.0 / (bpol(-nth:nth) * invRpol(-nth:nth)) if (debug) write(6,*) "eikcoefs: end gradients" da_bish = (1 / rmajor**2) + (bi / (bpol * rmajor**2))**2 db_bish = bi / (bpol * rmajor)**2 dc_bish = bi * (sin(th_bish) + rmajor * invRpol) / (-bpol * rmajor**4) da_bish = da_bish / trip db_bish = db_bish / trip dc_bish = dc_bish / trip call integrate(da_bish, outputs%theta, a_bish, ntgrid) call integrate(db_bish, outputs%theta, b_bish, ntgrid) call integrate(dc_bish, outputs%theta, c_bish, ntgrid) a_b = a_bish(nth) - a_bish(-nth) b_b = b_bish(nth) - b_bish(-nth) c_b = c_bish(nth) - c_bish(-nth) a_bish = - a_bish * bpol * rmajor b_bish = - b_bish * bpol * rmajor c_bish = - c_bish * bpol * rmajor ! find shat, pressure' and report them dp = geom%dpfun(pbar) ; di = geom%dbtori(pbar) if (abs(qval) > epsilon(0.)) then tmp = rhoc / (qval * twopi * drhodpsi) s_hat = tmp * (a_b * di + b_b * dp + c_b * 2) !< Calculate the shear from equib. if (local_eq) s_hat = surf%shat !Can't determine this by calculation in local_eq if (debug) write(6,*) "a_b (f'), b_b (p'), c_b terms =", a_b*di, b_b*dp, c_b*2 if (debug) write(6,*) 'qval, rho, drhodpsi, s_hat =', qval, rhoc, drhodpsi, s_hat else tmp = 1. ; s_hat = 0. end if write(fort11_unit,*) write(fort11_unit,*) 'Quantity, equilibrium value, value used' ! Decide how we want to set s_hat_new in the default case s_hat_new = s_hat select case (bishop) case (2, 3, 4, 5, 8, 9) !1, 6, 7 keep the initial value s_hat_new = surf%shat end select ! Decide how we want to set dp_new dp_new = dp select case (bishop) case (2) dp_new = p_prime_input case (3) pressure = geom%pfun(pbar) dp_new = -invLp_input * pressure * drhodpsi case (4, 6, 9) dp_new = 0.5 * beta_prime_input * drhodpsi case (5) dp_new = -alpha_input * drhodpsi / (qval**2 * rmaj) case (7) dp_new = dp * dp_mult case (8) dp_new = 0.5 * beta_prime_input * drhodpsi * dp_mult end select !Report on updated dp select case(bishop) case (2) write(fort11_unit,*) 'p_prime = ',dp,', ',dp_new if (writelots) write(*,*) 'p_prime = ',dp,', ',dp_new case (3) write(fort11_unit,*) '1/L_p = ',-dp / (drhodpsi * pressure), ', ', & -dp_new / (drhodpsi *pressure) if (writelots) write(*,*) '1/L_p = ',-dp / (drhodpsi * pressure), ', ', & -dp_new / (drhodpsi * pressure) case (5) write(fort11_unit,*) 'alpha = ', -qval**2 * rmaj * dp / drhodpsi,', ', & -qval**2 * rmaj * dp_new / drhodpsi if (writelots) write(*,*) 'alpha = ', -qval**2 * rmaj * dp / drhodpsi,', ', & -qval**2 * rmaj * dp_new / drhodpsi case default write(fort11_unit,*) 'd beta/d rho = ',2 * dp / drhodpsi,', ',2 * dp_new / drhodpsi if (writelots) write(*,*) 'd beta/d rho = ',2 * dp / drhodpsi,', ',2 * dp_new / drhodpsi end select write(fort11_unit,*) 's_hat ',s_hat,', ',s_hat_new if (writelots) write(*,*) 's_hat ',s_hat,', ',s_hat_new write(fort11_unit,*) if ((.not. local_eq) .and. (bishop == 1 .or. bishop > 9)) then di_new = di else di_new = (s_hat_new / tmp - 2 * c_b - dp_new * b_b) / a_b end if ! Calculate grad(S) ! here the coordinates are rho and l call gradl(ltheta, seik, dSdl, twopi * qval, nth) call periodic_copy(dSdl, 0., nperiod) grads(:,1) = a_bish * di_new + b_bish * dp_new + c_bish * 2 grads(:,2) = dSdl ! Calculate grad(Z) ! obtain Z(theta) and hence dZ/dl call Ztgrid(rgrid, outputs%theta, Zoftheta) call periodic_copy(Zoftheta, 0., nperiod) call gradl(ltheta, Zoftheta, dZdl, 0., nth) call periodic_copy(dZdl, 0., nperiod) ! gradztot(:,1) is grad Z . rhohat gradztot(:,1) = crpgrad(:,2) / outputs%grho ! gradztot(:,2) is grad Z . lhat gradztot(:,2) = dZdl ! gradztot(:,3) is grad Z . phihat gradztot(:,3) = 0.0 shat = s_hat_new ! Output dqdrp = s_hat_new * qval * drhodrp / rhoc ! Output dbetadrho = 2 * dp_new * dpsidrho ! Output call bishop_gradB(bmag, bpol, invRpol, th_bish, ltheta, rmajor, dp_new, bi, gradbtot) if (save_bishop) then call open_output_file(bishop_unit, '.bishop') write(bishop_unit, fmt="('#di, di_new = ',F12.7,' ',F12.7)") di, di_new write(bishop_unit, fmt="('#dp, dp_new = ',F12.7,' ',F12.7)") dp, dp_new write(bishop_unit, fmt="('#bi =', F12.7)") bi write(bishop_unit, fmt="('#drdpsi =', F12.7)") drhodpsi write(bishop_unit, fmt="('#a_b (fp)=', F20.10, ' ,b_b (pp)=', F20.10,' ,c_b =', F20.10,' ,s_hat =', F20.10)") a_b*di, b_b*dp, c_b, s_hat write(bishop_unit, fmt ="('# theta, a_bish, b_bish, c_bish, rmajor, Z, Bpol, invRpol, th_bish')") do i = -ntgrid, ntgrid write(bishop_unit, '( 10(2X, F25.17) )') outputs%theta(i), a_bish(i), b_bish(i), & c_bish(i), rmajor(i), Zoftheta(i), bpol(i), invRpol(i), th_bish(i) end do call close_output_file(bishop_unit) end if end subroutine bishop_grads !> @Note local_eq _must_ be false if we end in this routine. We make use of this. subroutine bishop_large_aspect_ratio_approximation(rhoc, dr, ntheta) real, intent(in) :: rhoc, dr integer, intent(in) :: ntheta real :: delrp, rp1, rp2, qval1, qval2, bi1, bi2, dqdrho, pbar_1, pbar_2 real, dimension(:), allocatable :: seik1, seik2, dsdrp integer :: itot, k ! rp derivative delrp = dr / drhodrp rp1 = rp - delrp ; rp2 = rp + delrp pbar_1 = pbarfun(rp1, outputs%theta(0)) ; pbar_2 = pbarfun(rp2, outputs%theta(0)) ! Calculate modified btori -- note we use the fact that rgrid == rp in ! seikon to avoid having to form rgrid1/rgrid2 here. bi1 = geom%btori(pbar_1) ; bi2 = geom%btori(pbar_2) allocate(seik1(-ntgrid:ntgrid), seik2(-ntgrid:ntgrid), dsdrp(-ntgrid:ntgrid)) call seikon(rp1, outputs%theta, invR, qval1, bi1, nperiod, seik1) call seikon(rp2, outputs%theta, invR, qval2, bi2, nperiod, seik2) dqdrho = (qval2 - qval1) / (2 * dr) shat = dqdrho * rhoc / qval ! Output dqdrp = dqdrho * drhodrp ! Output dbetadrho = 2 * geom%dpfun(pbar) / drhodrp ! Output if (writelots) then write (fort11_unit,*) 'dqdrho=',dqdrho,' qval1,2= ',qval1, qval2 write (fort11_unit,*) 's_hat= ',shat write (fort11_unit,*) 'i,rgrid,rgrid1,rgrid2:' do i = -nth, nth write (fort11_unit,*) i, rgrid(i), rp1, rp2 end do end if do i = -nth, nth dsdrp(i) = (seik2(i) - seik1(i)) / (2 * delrp) end do call periodic_copy(dsdrp, -twopi * dqdrp, nperiod) ! compute gradient of S ! this is really gradient of alpha ! note that coordinates are R and Z (not rho and l, as in bishop case) do k = -nperiod+1, nperiod-1 do i = -nth, nth itot = i + k * ntheta grads(itot, :) = rpgrad(i, :) * dsdrp(itot) + thgrad(i, :) * dsdthet(i) ! Output end do end do ! gradztot(:,1) is grad Z . Rhat ! Output gradztot(:,1) = 0. ! gradztot(:,2) is grad Z . Zhat ! Output gradztot(:,2) = 1. ! gradztot(:,3) is grad Z . phihat ! Output gradztot(:,3) = 0. call geom%gradient(rgrid, outputs%theta, gradbtot(:, 1:2), 'B', rp, nth, ntgrid, bishop /= 0) gradbtot(:, 3) = 0.0 end subroutine bishop_large_aspect_ratio_approximation !> FIXME : Add documentation subroutine alloc_local_arrays(n) implicit none integer, intent(in) :: n allocate(seik (-n:n), & dsdthet (-n:n), & invR (-n:n), & trip (-n:n), & rgrid (-n:n), & rmajor (-n:n)) allocate(thgrad (-n:n,2), & rpgrad (-n:n,2), & crpgrad (-n:n,2), & grads (-n:n,2)) allocate(gradrptot (-n:n,3), & gradstot (-n:n,3), & gradztot (-n:n,3), & gradthtot (-n:n,3), & ! MAB gradbtot (-n:n,3)) end subroutine alloc_local_arrays end subroutine eikcoefs !> Allocate the geom instance to the correct type pure subroutine allocate_geom(geom) use geq, only: geq_type use leq, only: leq_type use ceq, only: ceq_type use peq, only: peq_type, teq_type use eeq, only: eeq_type, gs2d_type use deq, only: deq_type use ideq, only: ideq_type use geo_utils, only: abstract_geo_type implicit none class(abstract_geo_type), allocatable, intent(in out) :: geom if (gen_eq) then allocate(geq_type :: geom) else if (chs_eq) then allocate(ceq_type :: geom) else if (gs2d_eq) then allocate(gs2d_type :: geom) else if (efit_eq) then allocate(eeq_type :: geom) else if (transp_eq) then allocate(teq_type :: geom) else if (ppl_eq) then allocate(peq_type :: geom) else if (dfit_eq) then allocate(deq_type :: geom) else if (idfit_eq) then allocate(ideq_type :: geom) else if (local_eq) then allocate(leq_type :: geom) end if end subroutine allocate_geom pure subroutine calculate_metric_terms(gradstot, gradrptot, gradthtot, dpsidrho, dqdrp, & qval, rhoc, bmag, force_sym, gds2, gds21, gds22, gds23, gds24, gds24_noq) real, dimension(-ntgrid:, :), intent(in) :: gradstot, gradrptot, gradthtot real, intent(in) :: dpsidrho, dqdrp, qval, rhoc real, dimension(-ntgrid:), intent(in) :: bmag logical, intent(in) :: force_sym real, dimension(-ntgrid:), intent(out) :: gds2, gds21, gds22, gds23, gds24, gds24_noq real, dimension(:), allocatable :: gdsdum1, gdsdum2, gdsdum3, gdsdum4, gdsdum5 allocate(gdsdum1(-ntgrid:ntgrid), gdsdum2(-ntgrid:ntgrid), gdsdum3(-ntgrid:ntgrid)) call dottgridf(gradstot, gradstot, gdsdum1) call dottgridf(gradstot, gradrptot, gdsdum2) call dottgridf(gradrptot, gradrptot, gdsdum3) gds2 = gdsdum1 *dpsidrho**2 gds21 = gdsdum2*dqdrp *dpsidrho**2 gds22 = gdsdum3*dqdrp*dqdrp*dpsidrho**2 ! v_E . grad theta = i*ky*phi*(rho_{r}*v_{t,r}/a**2)*(gds23 + gsd24*theta0) ! calculation of gds23 and gds24 assumes rp=rho, which is not necessarily true ! for numerical equilibria allocate(gdsdum4(-ntgrid:ntgrid), gdsdum5(-ntgrid:ntgrid)) call dottgridf(gradrptot, gradthtot, gdsdum4) !grad rho dot grad theta call dottgridf(gradstot, gradthtot, gdsdum5) !grad alpha dot grad theta gds23 = (gdsdum4 * gds2 - gdsdum5 * gdsdum2 * dpsidrho**2) / bmag**2 gds24_noq = (gdsdum4 * gdsdum2 - gdsdum5 * gdsdum3) * (dpsidrho / bmag)**2 gds24 = gds24_noq * dqdrp gds24_noq = gds24_noq * qval / rhoc if (force_sym) then call sym(gds2, 0, ntgrid) call sym(gds21, 1, ntgrid) call sym(gds22, 0, ntgrid) ! not sure about the following yet -- MAB ! call sym(gds23, 0, ntgrid); call sym(gds24, 1, ntgrid) end if end subroutine calculate_metric_terms subroutine calculate_surface_area(theta, gradpar, bmag, grho, dpsidrho, rhoc, nth, & surfarea, dvdrhon) use constants, only: twopi use file_utils, only: open_output_file, close_output_file real, dimension(-ntgrid:), intent(in) :: theta, gradpar, bmag, grho real, intent(in) :: dpsidrho, rhoc integer, intent(in) :: nth real, intent(out) :: surfarea, dvdrhon real, dimension(:), allocatable :: jacob, ds, dummy integer :: i, fort25_unit ! compute grad rho*Jacobian of the transformation to the new coordinates: ! Axisymmetry assumed to get this (nu indep. of phi): allocate(jacob(-ntgrid:ntgrid), ds(-ntgrid:ntgrid), dummy(-ntgrid:ntgrid)) do i = -nth, nth jacob(i) = abs(dpsidrho / (gradpar(i)*bmag(i))) ds(i) = grho(i) * jacob(i) end do ! compute surface area from A=Int(R sqrt(r**2+drdtheta**2) dtheta dphi) ! does not work with equal_arc grid ! compute the surface area from A=Int(J |grad rho| dtheta dalpha) call integrate(ds, theta, dummy, nth) surfarea = twopi * (dummy(nth) - dummy(-nth)) ! compute dV/drhon = 2 pi Int(J dtheta) call integrate(jacob, theta, dummy, nth) dvdrhon = twopi * (dummy(nth) - dummy(-nth)) call open_output_file(fort25_unit,".fort.25") write(fort25_unit,*) rhoc, f_trap(bmag(-nth:nth), theta(-nth:nth), jacob(-nth:nth)) call close_output_file(fort25_unit) end subroutine calculate_surface_area !> Computes something like x cross y dot grad zeta. pure subroutine tripprod2dtgrid(x, y, invR, val) implicit none real, dimension(-ntgrid:), intent (in) :: invR real, dimension(-ntgrid:), intent (out) :: val real, dimension(-ntgrid:, :), intent (in) :: x, y integer :: i ! factor of 1/R from grad zeta do i = -nth, nth val(i) = (x(i, 1) * y(i, 2) - x(i, 2) * y(i, 1)) * invR(i) end do end subroutine tripprod2dtgrid !> Calculate rgrid for the case of eeq and deq: in this case, rgrid !> is the distance to the magnetic axis as a function of theta. subroutine rtgrid(rgrid, rp, theta) implicit none real, dimension(-ntgrid:), intent (out) :: rgrid real, intent (in) :: rp real, dimension(-ntgrid:), intent (in) :: theta integer :: i do i = -nth, nth rgrid(i) = geom%rfun(rp, theta(i)) end do end subroutine rtgrid !> Calculates rmajor = R(theta), where R is the distance !> to the central axis (rmajor not to be confused with rmaj, r_geo etc) !> rgrid is the submodule radial grid (r_circ for EFIT & deq, rp otherwise) subroutine rmajortgrid(rgrid, theta, rmajor) implicit none real, dimension(-ntgrid:), intent(in) :: rgrid, theta real, dimension(-ntgrid:), intent(out) :: rmajor integer :: i do i = -nth, nth rmajor(i) = geom%Rpos(rgrid(i), theta(i)) end do end subroutine rmajortgrid !> Calculates invR = 1/R(theta), where R is the distance !> to the central axis (rmajor not to be confused with rmaj, r_geo etc) !> rgrid is the submodule radial grid (r_circ for EFIT & deq, rp otherwise) subroutine invRtgrid(rgrid, theta, invR) implicit none real, dimension(-ntgrid:), intent(in) :: rgrid, theta real, dimension(-ntgrid:), intent(out) :: invR integer :: i do i = -nth, nth invR(i) = geom%invR(rgrid(i), theta(i)) end do end subroutine invRtgrid !> Calculates the vertical height, Z, on the theta grid subroutine Ztgrid(rgrid, theta, Zoftheta) implicit none real, dimension(-ntgrid:), intent(in) :: rgrid, theta real, dimension(-ntgrid:), intent(out) :: Zoftheta integer :: i do i = -nth, nth Zoftheta(i) = geom%Zpos(rgrid(i), theta(i)) end do end subroutine Ztgrid !> FIXME : Add documentation pure subroutine bvectortgrid(invR, nth, gradrptot, dpsidrp, bi, bvector) integer, intent (in) :: nth real, dimension(-ntgrid: ), intent (in) :: invR real, dimension(-ntgrid:,:), intent (in) :: gradrptot real, dimension(-ntgrid:,:), intent (out) :: bvector real, intent (in) :: dpsidrp, bi integer :: i do i = -nth, nth bvector(i,1) = -dpsidrp * gradrptot(i,2) * invR(i) bvector(i,2) = dpsidrp * gradrptot(i,1) * invR(i) bvector(i,3) = bi * invR(i) end do end subroutine bvectortgrid !> FIXME : Add documentation pure subroutine gradstottgrid (invR, grads, gradstot) real, dimension(-ntgrid:), intent (in) :: invR real, dimension(-ntgrid:, :), intent (in) :: grads real, dimension(-ntgrid:, :), intent (out) :: gradstot gradstot(:, 1:2) = grads if (use_large_aspect) then gradstot(:, 3) = 0. else gradstot(:, 3) = invR end if end subroutine gradstottgrid !> Calculate the cross product c = a x b on the `[-nth, nth]` domain pure subroutine crosstgrid(a, b, c) real, dimension(-ntgrid:, :), intent (in) :: a, b real, dimension(-ntgrid:, :), intent (out) :: c integer :: i do i = -nth, nth c(i, 1) = a(i, 2) * b(i, 3) - b(i, 2) * a(i, 3) c(i, 2) = a(i, 3) * b(i, 1) - b(i, 3) * a(i, 1) c(i, 3) = a(i, 1) * b(i, 2) - b(i, 1) * a(i, 2) end do end subroutine crosstgrid !> Calculate the dot product c = a.b where a, b and c are on the `[-ntgrid, ntgrid]` domain pure subroutine dottgridf(a, b, c) real, dimension(-ntgrid:,:), intent(in) :: a, b real, dimension(-ntgrid:), intent(out) :: c c = a(:, 1) * b(:, 1) + a(:, 2) * b(:, 2) + a(:, 3) * b(:, 3) end subroutine dottgridf !> Returns the diameter of the surface labelled by rp. real function diameter(rp) implicit none real, intent(in) :: rp diameter = geom%diameter(rp) !Used to allow rp<rpmin here by just returning zero end function diameter !> FIXME : Add documentation !> Only used with ~undocumented `irho = 4` option. Ideally would update to take pbar real function rhofun (rp) implicit none real, intent (in) :: rp real :: pbar pbar = psifun(rp) rhofun = geom%rhofun(pbar) end function rhofun !> Calculates the normalised radial coordinate pure real function psifun (rp) implicit none real, intent (in) :: rp psifun = min(1.0, max(0.0, (rp - rpmin) / (rpmax - rpmin))) end function psifun real function pbarfun(r, theta) implicit none real, intent(in) :: r, theta pbarfun = psifun(geom%psi(r, theta)) end function pbarfun !> FIXME : Add documentation real function pbarofrho(rho) implicit none real, intent (in) :: rho pbarofrho = psifun(rpofrho(rho)) end function pbarofrho !> FIXME : Add documentation subroutine drift(rgrid, rp, gradstot, gradrptot, gradztot, gradbtot, dqdrp, dpsidrho, & dpsidrp, bmag, invR, theta, dbetadrho, bi, force_sym, nperiod, & gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0) implicit none real, dimension (-ntgrid:), intent (in) :: rgrid, bmag, invR, theta real, dimension (-ntgrid:, :), intent (in) :: gradstot, gradrptot, gradztot, gradbtot real, intent (in) :: rp, dqdrp, dpsidrho, dpsidrp, dbetadrho, bi logical, intent(in) :: force_sym integer, intent(in) :: nperiod real, dimension(-ntgrid:ntgrid), intent(out) :: gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0 real, dimension (-ntgrid:ntgrid, 3) :: pgradtot, dummy, curve, bvector real, dimension (-ntgrid:ntgrid) :: gbdrift1, gbdrift2, cvdrift1, cvdrift2, cdrift1, cdrift2 real :: dp if (use_large_aspect) then gbdrift = (-sin(theta) * gradstot(:,1) - cos(theta) * gradstot(:,2)) * invR gbdrift = 2 * dpsidrho * gbdrift gbdrift0 = (-sin(theta) * gradrptot(:,1) - cos(theta) * gradrptot(:,2)) * invR gbdrift0 = 2 * dpsidrho * gbdrift0 * dqdrp cvdrift = gbdrift ; cvdrift0 = gbdrift0 else ! compute magnetic field call bvectortgrid(invR, nth, gradrptot, dpsidrp, bi, bvector) dp = dbetadrho / (2 * dpsidrho) call crosstgrid(bvector, gradbtot, dummy) !B cross Grad B call periodic_copy(dummy(:, 1), 0.0, nperiod) call periodic_copy(dummy(:, 2), 0.0, nperiod) call periodic_copy(dummy(:, 3), 0.0, nperiod) call dottgridf(dummy, gradstot, gbdrift1) call dottgridf(dummy, gradrptot, gbdrift2) ! create curvature if (bishop == 0) then call geom%gradient(rgrid, theta, pgradtot(:, 1:2), 'R', rp, nth, ntgrid, bishop /= 0) pgradtot(-nth:nth, 3) = 0.0 else pgradtot = gradrptot * dpsidrp * dp !Grad rho * dp/drho end if curve = gradbtot curve(-nth:nth, 1) = curve(-nth:nth, 1) + pgradtot(-nth:nth, 1) / bmag(-nth:nth) curve(-nth:nth, 2) = curve(-nth:nth, 2) + pgradtot(-nth:nth, 2) / bmag(-nth:nth) call crosstgrid(bvector, curve, dummy) !B cross [Grad B + Grad P / B] call periodic_copy(dummy(:, 1), 0.0, nperiod) call periodic_copy(dummy(:, 2), 0.0, nperiod) call periodic_copy(dummy(:, 3), 0.0, nperiod) call dottgridf(dummy, gradstot, cvdrift1) call dottgridf(dummy, gradrptot, cvdrift2) ! factor of 2 appears below because will later multiply ! by wunits, which has a factor of 1/2 built-in gbdrift = 2 * dpsidrho * gbdrift1 / bmag**3 gbdrift0 = 2 * dpsidrho * gbdrift2 * dqdrp / bmag**3 cvdrift = 2 * dpsidrho * cvdrift1 / bmag**3 cvdrift0 = 2 * dpsidrho * cvdrift2 * dqdrp / bmag**3 end if ! coriolis drift term is -2*vpa*om_phi/Omega * (bhat x (zhat x bhat)) . grad (g+q<phi>F_0/T) ! here we calculate 2*(zhat . grad (alpha+q*theta0))/(B/Bref) ! this is the Zhat . grad(alpha) part of Zhat . grad(S) ! note that a factor of k_y = (n0/a) (drho/dpsiN) has been factored out call dottgridf(gradztot, gradstot, cdrift1) ! this is the Zhat . grad(q) part of Zhat . grad(S) ! note that a factor of k_y * theta0 has been factored out call dottgridf(gradztot, gradrptot, cdrift2) ! factor of 4 (instead of 2) below because will later mutliply ! by wunits, which has a built-in factor of 1/2 ! cdrift0 is part of coriolis drift dependent on theta0 cdrift = -4 * dpsidrho * cdrift1 / bmag cdrift0 = -4 * dpsidrho * cdrift2 * dqdrp / bmag if (force_sym) then call sym(gbdrift, 0, ntgrid) call sym(cvdrift, 0, ntgrid) call sym(gbdrift0, 1, ntgrid) call sym(cvdrift0, 1, ntgrid) ! should maybe add in cdrift and cdrift0 here -- MAB end if end subroutine drift !> Forces symmetry in theta in the passed array. If isign == 0 then !> `a(+/- theta)` is replaced by the average of `a(+/-theta)`, otherwise !> `a(+/- theta)` is replaced by `+/- (a(theta)-a(-theta))`. pure subroutine sym(a, isign, ntgrid) implicit none integer, intent(in) :: isign, ntgrid real, dimension(-ntgrid:), intent (in out) :: a integer :: i if (isign == 0) then do i = 1, ntgrid a(i) = 0.5 * (a(i) + a(-i)) a(-i) = a(i) end do else do i = 1, ntgrid a(i) = 0.5 * (a(i) - a(-i)) a(i) = -a(-i) end do a(0) = 0. end if end subroutine sym !> FIXME: Add documentation. Only used for bishop == 0 testing subroutine seikon(rp, theta, invR, qval, bi, nperiod, seik) real, intent(in) :: rp real, dimension(-ntgrid:), intent(in) :: theta, invR real, intent(in out) :: qval real, intent(in) :: bi integer, intent(in) :: nperiod real, dimension(-ntgrid:), intent(out) :: seik real :: dpsidrp real, dimension(-ntgrid:ntgrid, 2) :: thgrad, rpgrad real, dimension(-ntgrid:ntgrid) :: trip, dsdthet, rgrid rgrid = rp !<This is not right for cart types, but ok for those compatible with bishop=0. call geom%gradient(rgrid, theta, rpgrad, 'P', rp, nth, ntgrid, bishop /= 0) call geom%gradient(rgrid, theta, thgrad, 'T', rp, nth, ntgrid, bishop /= 0) call tripprod2dtgrid(rpgrad, thgrad, invR, trip) ; call periodic_copy(trip, 0., nperiod) call eikonal(theta, invR, trip, qval, bi, seik, dsdthet, dpsidrp) end subroutine seikon !> FIXME : Add documentation pure subroutine eikonal(theta, invR, trip, qval, bi, seik, dsdthet, dpsidrp) use constants, only: twopi implicit none real, dimension(-ntgrid:), intent(in) :: theta, invR, trip real, dimension(-ntgrid:), intent(out) :: seik, dsdthet real, intent(in out) :: qval real, intent(in) :: bi real, intent(out) :: dpsidrp dsdthet = -bi * invR**2 / trip call integrate(dsdthet, theta, seik, ntgrid) if (local_eq) then dpsidrp = -(seik(nth) - seik(-nth)) / (twopi * qval) seik = seik / dpsidrp dsdthet = dsdthet / dpsidrp else dpsidrp = 1. qval = -(seik(nth) - seik(-nth)) / (twopi) end if end subroutine eikonal !> FIXME : Add documentation real function rpofrho(rho) implicit none real, intent(in) :: rho real, parameter :: xerrbi = 1.e-4, xerrsec = 1.e-8 interface real function f(x_) implicit none real, intent(in) :: x_ end function f end interface procedure(f), pointer :: func real :: a, b, fval integer :: ier ier = 123456 ! Default to indicating error a = rpmin ; b = rpmax if (irho == 1) then func => phi ; fval = rho * rho * func(b) else if (irho == 2) then func => diameter ; fval = rho * func(b) else if (irho == 3) then func => psifun ; fval = rho * func(b) else if (irho == 4) then func => rhofun ; fval = rho * func(b) end if call root(func, fval, a, b, xerrbi, xerrsec, ier, rpofrho) if (debug) then write (*,*) "Values in root: " write (*,*) 'rho: ', rho write (*,*) 'fval: ', fval write (*,*) 'rpofrho: ', rpofrho write (*,*) 'rpmin: ', rpmin write (*,*) 'rpmax: ', rpmax end if if (ier > 0) write(*, *) 'error in rpofrho, rho=', rho, rpofrho, ier end function rpofrho !> Integrates arg(grid) from i = 0 (generally theta = 0) to either !> ends of the grid (usually +/- pi) storing the cumulative integral !> in the output `ans`. pure subroutine integrate(arg, grid, ans, n) real, dimension(-ntgrid:), intent (in) :: arg, grid real, dimension(-ntgrid:), intent (out) :: ans integer, intent(in) :: n integer :: i ans = 0. do i = 1, n ans(i) = ans(i-1) + 0.5 * (grid(i) - grid(i-1)) * (arg(i) + arg(i-1)) end do do i = -1, -n, -1 ans(i) = ans(i+1) + 0.5 * (grid(i) - grid(i+1)) * (arg(i) + arg(i+1)) end do end subroutine integrate !> Find soln where f(soln) ~= fval. Works best if a and b bracket the root. !> We can then first bisect the bracket to narrow in on the solution. !> Following this we use the secant method to find the solution. !> The input xerrbi influences the number of bisection iterations whilst !> xerrsec is used to set the stopping tolerance on the secant method. subroutine root(f, fval, a, b, xerrbi, xerrsec, ier, soln) implicit none interface real function f(x_) implicit none real, intent(in) :: x_ end function f end interface real, intent(in) :: fval, a, b, xerrbi, xerrsec real, intent(out) :: soln integer, intent(out) :: ier real :: a1, b1, f1, f2, f3, trial, aold integer :: i, niter logical :: skip_bisection ier = 0 ; a1 = a ; b1 = b f1 = f(a1) - fval ; f2 = f(b1) - fval skip_bisection = xerrbi <= 0 !Could be xerrbi <= abs(b1 - a1)? if (.not. skip_bisection .and. f1 * f2 > 0.) then write(*, *) 'f1 and f2 have same sign in bisection routine root' write(*, *) 'a1,f1,b1,f2=', a1, f1, b1, f2 write(*, *) 'fval=', fval ier = 1 skip_bisection = .true. end if if (.not. skip_bisection) then ! Presumably we require |b1 - a1| >= xerrbi here to avoid niter < 1? niter = 1 + int(log(abs(b1 - a1) / xerrbi) / log(2.)) ! Repeatedly bisect between limits, swapping one limit ! for the bisection point each iteration. Helps to narrow ! down the bracketing region substantially. do i = 1, niter trial = 0.5 * (a1 + b1) f3 = f(trial) - fval if (f3 * f1 > 0.) then a1 = trial f1 = f3 else b1 = trial f2 = f3 end if end do end if !Swap so |f1| < |f2| if( abs(f1) > abs(f2) ) then f3 = f1 ; f1 = f2 ; f2 = f3 aold = a1 ; a1 = b1 ; b1 = aold end if ! start secant method do i = 1, 10 aold = a1 ; f3 = f2 - f1 if (abs(f3) < 1.e-11) f3 = 1.e-11 !Should we exit here? Saying f2~f1 a1 = a1 - f1 * (b1 - a1) / f3 !Update guess for root. f2 = f1 ; b1 = aold ; f1 = f(a1) - fval if (abs(a1 - b1) < xerrsec) exit !Root found closely enough so exit end do soln = a1 ! Detect instances where we have not converged to a solution if (abs(f1) > xerrsec) ier = ier + 2 end subroutine root !> FIXME : Add documentation real function phi(rp) implicit none real, intent (in) :: rp integer, parameter :: nimax = 200 real :: pbar, pb(nimax), dpb, qi, qnext integer :: i pbar = psifun(rp) ! Construct a uniformly space grid from 0 to pbar dpb = pbar / float(nimax - 1) do i = 1, nimax pb(i) = dpb * float(i - 1) end do ! Integrate area under q(pb) curve from pb = 0 to pbar phi = 0. ; qi = geom%qfun(pb(1)) do i = 1, nimax - 1 qnext = geom%qfun(pb(i+1)) phi = phi + 0.5 * dpb * (qnext + qi) qi = qnext end do phi = phi * (rpmax - rpmin) end function phi !> Given an initial theta and gradpar (=dtheta/dl) compute a new !> theta grid, returned in arcl, which has gradpar = dg/dl = constant. !> @note The input gradpar is modified on output. pure subroutine arclength (ntheta, nperiod, theta, gpar, arcl) use constants, only: twopi, pi implicit none integer, intent (in) :: ntheta, nperiod real, dimension(-ntgrid:), intent (out) :: arcl real, dimension(-ntgrid:), intent (in) :: theta real, dimension(-ntgrid:), intent (in out) :: gpar integer :: nth, j real :: arcfac nth = ntheta / 2 arcl(-nth) = 0. do j = -nth, nth-1 arcl(j+1) = arcl(j) + 0.5 * (theta(j+1) - theta(j)) * (1. / gpar(j) + 1. / gpar(j+1)) end do arcfac = 1. / arcl(nth) do j = -nth, nth arcl(j) = arcl(j) * twopi * arcfac - pi end do call periodic_copy(arcl, twopi, nperiod) gpar = twopi * arcfac end subroutine arclength !> FIXME : Add documentation subroutine loftheta(rgrid, theta, ltheta) implicit none real, dimension (-ntgrid:), intent(in) :: rgrid, theta real, dimension (-ntgrid:), intent(out) :: ltheta integer :: i real :: rold, told, rpos_old, zpos_old, rpos_new, zpos_new real :: r, thet, rpos_orig, zpos_orig ltheta(0) = 0. rold = rgrid(0) ; told = theta(0) rpos_orig = geom%Rpos(rold, told) ; zpos_orig = geom%Zpos(rold, told) rpos_new = rpos_orig ; zpos_new = zpos_orig do i = 1, nth rpos_old = rpos_new ; zpos_old = zpos_new r = rgrid(i) ; thet = theta(i) rpos_new = geom%Rpos(r, thet) ; zpos_new = geom%Zpos(r, thet) ltheta(i) = ltheta(i-1) + hypot((Rpos_new - Rpos_old), (Zpos_new - Zpos_old)) end do rpos_new = rpos_orig ; zpos_new = zpos_orig do i = -1, -nth, -1 rpos_old = rpos_new ; zpos_old = zpos_new r = rgrid(i) ; thet = theta(i) rpos_new = geom%Rpos(r, thet) ; zpos_new = geom%Zpos(r, thet) ltheta(i) = ltheta(i+1) - hypot((Rpos_new - Rpos_old), (Zpos_new - Zpos_old)) end do end subroutine loftheta !> Compute df / dltheta using a periodic central difference. The !> input `ext` can be used to control the effective jump in f !> between both ends of the "periodic" grid. pure subroutine gradl (ltheta, f, dfdl, ext, n) implicit none real, dimension (-ntgrid:), intent(in) :: ltheta, f real, dimension (-ntgrid:), intent(out) :: dfdl real, intent(in) :: ext integer, intent(in) :: n integer :: i dfdl(-n) = (f(n-1) + ext - f(-n+1)) / & (ltheta(n-1) - ltheta(-n+1) - ltheta(n) + ltheta(-n)) do i = -n+1, n-1 dfdl(i) = (f(i+1) - f(i-1)) / (ltheta(i+1) - ltheta(i-1)) end do dfdl(n) = dfdl(-n) end subroutine gradl !> Finds the angle between the surface tangent vector in the poloidal plane !> and the direction of increasing major radius. pure subroutine th_bishop(rpgrad, th_bish, nth) use constants, only: twopi implicit none integer, intent (in) :: nth real, dimension (-ntgrid:, :), intent (in) :: rpgrad real, dimension (-ntgrid:), intent (out) :: th_bish real, dimension (-nth:nth) :: magrp real, dimension (-nth:nth, 2) :: tvec integer :: i real :: tvec_dot_gradR magrp = hypot(rpgrad(-nth:nth, 1), rpgrad(-nth:nth, 2)) ! tvec is the tangent vector in the surface of the poloidal plane tvec(:, 1) = rpgrad(-nth:nth, 2) / magrp tvec(:, 2) = -rpgrad(-nth:nth, 1) / magrp ! need to find grad(R) and dot it with tvec. but grad(R) is ! simple -- it is just (1,0) so simply return acos(tvec(:,1)) with ! corrections for left-hand quadrants do i = -nth, nth tvec_dot_gradR = min(1., max(-1., tvec(i, 1))) if (tvec(i, 2) > 0.) then th_bish(i) = twopi - acos(tvec_dot_gradR) else th_bish(i) = acos(tvec_dot_gradR) end if end do end subroutine th_bishop !> FIXME : Add documentation pure subroutine B_pol(invR, rpgrad, dpsidrp, bpol, nth) implicit none integer, intent (in) :: nth real, intent(in) :: dpsidrp real, dimension (-ntgrid:), intent (in) :: invR real, dimension (-ntgrid:, :), intent (in) :: rpgrad real, dimension (-ntgrid:) , intent (out) :: bpol bpol(-nth:nth) = hypot(rpgrad(-nth:nth, 1), rpgrad(-nth:nth, 2)) * invR(-nth:nth) * dpsidrp end subroutine B_pol !> FIXME : Add documentation pure subroutine B_mod(invR, bpol, bmag, bi, nth) implicit none integer, intent (in) :: nth real, dimension (-ntgrid:), intent (in) :: invR, bpol real, dimension(-ntgrid:), intent (out) :: bmag real, intent(in) :: bi bmag(-nth:nth) = hypot((bi * invR(-nth:nth)), bpol(-nth:nth)) end subroutine B_mod !> General calculation of the poloidal magnetic field radius of curvature !> (see section 2.3 of "GS2 analytic geometry specification") ! JRB !> @Note we actually calculate 1/Rpol as this is what we actually use. pure subroutine inv_Rpol(theta, th_bish, ltheta, invRpol, nth) use constants, only: twopi, pi implicit none integer, intent (in) :: nth real, dimension(-ntgrid:), intent (in) :: th_bish, theta, ltheta real, dimension(-ntgrid:), intent (out) :: invRpol real, dimension(-ntgrid:ntgrid) :: dthdl integer :: i real :: dth_bish, dth ! Calculate du/dtheta (where u is defined according to Miller PoP 5 973 (1998)) ! Note that we must properly account for discontinuities in u. u is th_bish do i = -nth + 1, nth - 1 dth_bish = th_bish(i + 1) - th_bish(i - 1) ; dth = theta(i + 1) - theta(i - 1) if ( dth_bish > pi) then invRpol(i) = (dth_bish - twopi) / dth else if(dth_bish < -pi) then invRpol(i) = (dth_bish + twopi) / dth else invRpol(i) = dth_bish / dth end if end do invRpol(nth) = (th_bish(-nth+1) - th_bish(nth-1)) / (theta(-nth+1) - theta(nth-1) + twopi) invRpol(-nth) = invRpol(nth) ! Calculate dtheta/dl (where l is the poloidal arc length) call gradl(ltheta, theta, dthdl, -twopi, nth) ! Scale by dthdl to complete calculation of 1/Rpol ! (the inverse radius of curvature of the poloidal magnetic field) invRpol(-nth:nth) = invRpol(-nth:nth) * dthdl(-nth:nth) end subroutine inv_Rpol !> FIXME : Add documentation subroutine finish_geometry implicit none if (allocated(geom)) then call geom%finalise deallocate(geom) end if end subroutine finish_geometry !> Initialse [[geometry:theta]] to a uniform grid with `2*((2*nperiod - 1)*(nt/2)) + 1` !> points. Also set [[geometry:nth]], [[geometry:ntheta]] and [[geometry::ntgrid]] !> consistently; importantly [[geometry:nth]] is forced to be even. If `ntheta` is odd, then !> `nth == ntheta - 1`. !> !> `ntheta` is set equal to [[geometry:ntheta]] on exit subroutine init_uniform_theta_grid(ntheta, theta, nperiod, ntgrid, nth) use constants, only: pi implicit none integer, intent(in out) :: ntheta real, dimension(:), allocatable, intent(in out) :: theta integer, intent(in) :: nperiod integer, intent(out) :: ntgrid, nth integer :: i nth = ntheta / 2 ! Make sure ntheta is even, update argument and set ntgrid ntheta = nth * 2 ; ntgrid = (2 * nperiod - 1) * nth if (verb>3) write(6,*) "init_theta: allocated(theta),ntgrid,ntheta,nperiod=",& allocated(theta), ntgrid, ntheta, nperiod if (allocated(theta)) deallocate (theta) allocate(theta(-ntgrid:ntgrid)) theta = [ (i * pi / real(nth), i=-ntgrid, ntgrid) ] end subroutine init_uniform_theta_grid !> Copy the `[-nth,nth] = [-pi,pi]` section of the passed array into !> the `theta > +/- pi` parts, adding on `k * ext` (with k an integer !> labelling the sections/periods). Here we expect `ext == array(nth) - array(-nth)`. pure subroutine periodic_copy(a, ext, nperiod) implicit none real, dimension(-ntgrid:), intent(in out) :: a real, intent(in) :: ext integer, intent(in) :: nperiod integer :: i, k, itot, ntheta if (nperiod <= 1) return ! Note we end up enforcing a(nth)-a(-nth) = ext, but we should probably check this is ! also true on input to the routine ntheta = 2 * nth do k = -nperiod + 1, -1 do i = -nth, nth itot = i + k * ntheta a(itot) = a(i) + k * ext end do end do do k = 1, nperiod - 1 do i = -nth, nth itot = i + k * ntheta a(itot) = a(i) + k * ext end do end do end subroutine periodic_copy !> FIXME : Add documentation pure subroutine bishop_gradB(bmag, Bpol, invRpol, th_bish, ltheta, rmajor, dp, bi, gradB) implicit none real, dimension(-ntgrid:), intent (in) :: bmag, bpol, invRpol real, dimension(-ntgrid:), intent (in) :: th_bish, ltheta, rmajor real, intent(in) :: dp, bi real, dimension(-ntgrid:,:), intent(out) :: gradB real :: bp integer :: i do i = -nth, nth bp = -bpol(i) gradB(i, 1) = (bp / bmag(i)) * (bp * invRpol(i) + dp * rmajor(i) & - bi**2 * sin(th_bish(i)) / (rmajor(i)**3 * bp)) !dBdrho end do call gradl(ltheta, bmag, gradB(:, 2), 0., nth) !dBdl gradB(:, 3) = 0.0 end subroutine bishop_gradB !> Perform minor checks on consistency of flags subroutine check(geq, eeq, peq, leq, deq, ideq, ceq, gs2deq, teq, bishop) use mp, only: mp_abort implicit none logical, intent(in) :: geq, eeq, peq, leq, deq, ideq, ceq, gs2deq, teq integer, intent(in) :: bishop integer :: nflag_set if (.not.(geq .or. ideq) .and. bishop == 0) & !Possibly also need use_large_aspect if bishop == 0? call mp_abort('bishop = 0 only compatible with geq or ideq',.true.) nflag_set = count([geq, ceq, eeq, gs2deq, peq, teq, leq, deq, ideq]) if (nflag_set == 0) then call mp_abort('Must set at least one equilibrium flag to true', .true.) else if (nflag_set > 1) then call mp_abort('More than one equilibrium flag to true, not allowed', .true.) end if end subroutine check !> Calculate the derivative of the flux label rho with respect to !> the internal flux label rp (the normalised poloidal flux) real function drho_drp(rp, dr) implicit none real, intent(in) :: rp, dr real :: rp1, rp2, rho1, rho2, phirpmax rp1 = rp * (1 - dr) ; rp2 = rp * (1 + dr) if (irho == 1) then phirpmax = phi(rpmax) rho1 = sqrt(abs(phi(rp1) / phirpmax)) ; rho2 = sqrt(abs(phi(rp2) / phirpmax)) else if (irho == 2) then rho1 = 0.5 * diameter(rp1) ; rho2 = 0.5 * diameter(rp2) else if (irho == 3) then rho1 = psifun(rp1) ; rho2 = psifun(rp2) else if (irho == 4) then rho1 = rhofun(rp1) ; rho2 = rhofun(rp2) end if drho_drp = (rho2 - rho1) / (rp2 - rp1) end function drho_drp !> FIXME : Add documentation | Only used for writelots output and if irho=1 real function drho_drhod(rp, drhodrp, dr) implicit none real, intent (in) :: drhodrp, rp, dr real :: rp1, rp2, rhod1, rhod2 ! compute d rho / d rb rp1 = rp * (1 - dr) ; rp2 = rp * (1 + dr) rhod1 = 0.5 * diameter(rp1) ; rhod2 = 0.5 * diameter(rp2) drho_drhod = drhodrp / ((rhod2 - rhod1) / (rp2 - rp1)) end function drho_drhod !> FIXME : Add documentation pure real function fluxavg(f, theta, jacob) implicit none real, dimension(-nth:nth), intent(in) :: f, theta, jacob real, dimension(-nth:nth) :: delth delth(-nth+1:nth-1) = 0.5 * (theta(-nth+2:nth) - theta(-nth:nth-2)) delth(-nth) = 0.5 * (theta(-nth+1) - theta(-nth)) delth(nth) = 0.5 * (theta(nth) - theta(nth-1)) ! We might want to write this to use trapezoidal_integration instead ! of this sum form. Currently only used in [[f_trap]], which itself ! is only used here for one piece of information written to screen. ! As it doesn't get into the simulated system we probably don't mind ! too much beyond readability concerns. fluxavg = sum(f * jacob * delth) / sum(jacob * delth) end function fluxavg !> FIXME : Add documentation pure real function f_trap(b_mag, theta, jacob) implicit none real, dimension(-nth:nth), intent (in) :: b_mag, theta, jacob real :: ftu, ftl, havg, h2avg, h(-nth:nth) ! use expressions from Y. R. Lin-Liu and Bob Miller, coding from Q. P. Liu h = min(b_mag / maxval(b_mag),1.0) havg = fluxavg(h, theta, jacob) h2avg = fluxavg(h**2, theta, jacob) ftu = 1.0 - (h2avg / havg**2) * (1.0 - sqrt(1.0 - havg) * (1.0 + 0.5 * havg)) ftl = 1.0 - h2avg * fluxavg((1.0 - sqrt(1.0 - h) * (1.0 + 0.5 * h)) / h**2, theta, jacob) f_trap = 0.75 * ftu + 0.25 * ftl end function f_trap !> FIXME : Add documentation subroutine plotdata (rgrid, seik, grads, dpsidrho, dr, theta, & rplot, zplot, aplot, rprime, zprime, aprime) implicit none real, dimension (-ntgrid:), intent(in) :: rgrid, seik, theta real, dimension (-ntgrid:, :), intent(in) :: grads real, dimension (-ntgrid:), intent(out) :: rplot, zplot, aplot, rprime, zprime, aprime real, intent (in) :: dpsidrho, dr real :: rplus, rminus, dr_tot integer :: i ! At the present flux surface, calculate these three functions ! as a function of theta, together with their derivates with ! respect to rho: ! ! R, Z, alpha-phi ! ! Call these variables ! ! Rplot, Zplot, aplot ! ! and their derivatives: ! ! Rprime, Zprime, aprime ! ! BD: bug?? ! Is this correct for EFIT data? Need to check do i = -ntgrid, ntgrid Rplot(i) = geom%Rpos(rgrid(i), theta(i)) Zplot(i) = geom%Zpos(rgrid(i), theta(i)) end do aplot = seik do i = -ntgrid, ntgrid rplus = rgrid(i) * (1 + dr) ; rminus = rgrid(i) * (1 - dr) dr_tot = 2 * dr * rgrid(i) Rprime(i) = (geom%Rpos(rplus, theta(i)) - geom%Rpos(rminus, theta(i))) / dr_tot Zprime(i) = (geom%Zpos(rplus, theta(i)) - geom%Zpos(rminus, theta(i))) / dr_tot end do aprime = grads(:, 1) * dpsidrho end subroutine plotdata end module geometry