!> This module is a submodule of geometry which handles reading from the !> ascii EQDSK format output by EFIT, but also now by other codes. This !> file contains psi on an R,Z grid, as well as other quantities such as !> q (safety factor), I (a.k.a. f) and p (pressure) on a psi grid. !> !> The normalising field is set to the field on axis; the normalising length !> is set to the half-diameter of the LCFS. module eeq use geo_utils, only: abstract_eqfile_cart_geo_type, geo_input_type, sort implicit none private public :: eeq_type, gs2d_type type, abstract, extends(abstract_eqfile_cart_geo_type) :: eeq_base_type private contains procedure :: setup_special_splines, delete_special_splines procedure :: rcenter, dealloc_special_arrays end type eeq_base_type type, extends(eeq_base_type) :: eeq_type contains procedure :: read_and_set => efit_read_and_set end type eeq_type type, extends(eeq_base_type) :: gs2d_type contains procedure :: read_and_set => gs2d_read_and_set end type gs2d_type contains !> This subroutine reads a GS2D output file containing !> the axisymmetric magnetic field geometry on a rectangular !> domain defined by the coordinates (R,Z). subroutine gs2d_read_and_set(self, inputs) use constants, only: twopi use gs2d, only: read_gs2d, psi, b0, f, ippsi, nr, nz, p, ps use gs2d, only: psip, psmin, q, r0, rgrid, rmag, rsep, zgrid, zmag, zsep use mp, only: mp_abort implicit none class(gs2d_type), intent(in out) :: self type(geo_input_type), intent(in) :: inputs real :: rwid, rleft, zhei, rcentr real, dimension(:), allocatable :: spsi_bar real, dimension(:, :), allocatable :: sefit_psi integer :: nw_in, nh_in logical, parameter :: debug =.false. self%type_name = 'GS2DEQ' self%filename = trim(adjustl(inputs%eqfile)) if (debug) write(6,*) "gs2din: filename=",self%filename ! read GS2D datafile call read_gs2d(self%filename) if (debug) write(6,*) "gs2din: read_gs2d done, " if (debug) write(6,fmt='(T2,"psmin",T15,"psedge (Wb)",T27,"b0 (T)",T39,"ip(A)"/1p,4e16.8)') psmin,psip,b0,ippsi nw_in = nr ; nh_in = nz if (size(ps) /= nw_in) then write(6,*) 'gs2din: size(ps) (',size(ps),') /= nw (',nw_in,')' write(6,*) 'gs2din: => should fix the GS2D output file and try again' call mp_abort('gs2din: => should fix the GS2D output file and try again') end if allocate(spsi_bar(nw_in), sefit_psi(nw_in, nh_in)) rwid = maxval(rgrid) - minval(rgrid) ; zhei = maxval(zgrid) - minval(zgrid) rcentr = r0 ; rleft = minval(rgrid) self%R_mag = rmag ; self%Z_mag = zmag self%psi_0 = psmin ; self%psi_a = psip self%psi_0 = self%psi_0 / twopi ; self%psi_a = self%psi_a / twopi if (debug) write(6,*) "gs2din: psi_0, psi_a=", self%psi_0, self%psi_a ! divide GS2D psi by 2pi to get poloidal flux in (Wb/rad) sefit_psi = psi / twopi ! pbar is defined by pbar = (psi-psi_0)/(psi_a-psi_0). fp and q are ! functions of pbar spsi_bar = (ps - minval(ps)) / (maxval(ps) - minval(ps)) call self%shared_setup(nw_in, nh_in, spsi_bar, rgrid, zgrid, f, p, q, sefit_psi, & size(rsep), rsep, zsep, b0, rleft, rwid, minval(zgrid), zhei, inputs%big, .true.) end subroutine gs2d_read_and_set !> This subroutine reads an EFIT output file containing !> the axisymmetric magnetic field geometry on a rectangular !> domain defined by the coordinates (R,Z). It reads and stores !> the following quantities (among others). !> !> - eq_R: R grid !> - eq_Z: Z grid !> - fp: F on psibar grid !> - eqpsi_2d: psi on (R,Z) grid !> - R_mag: major radius of the magnetic axis !> - Z_mag: elevation of the magnetic axis !> - rwid: total width of the domain !> - rleft: position of leftmost point of domain !> - zhei: total height of domain subroutine efit_read_and_set(self, inputs) implicit none class(eeq_type), intent(in out) :: self type(geo_input_type), intent(in) :: inputs real :: xdum, rwid, rleft, zhei, rcentr, bcentr, zmid real, dimension(:), allocatable :: dummy, spsi_bar, sefit_R, sefit_Z, rbbbs, zbbbs, f, p, q real, dimension(:, :), allocatable :: sefit_psi character(10) :: char integer :: i, j, ndum, nw_in, nh_in, funit, nbbbs self%type_name = 'EFIT' self%filename = trim(adjustl(inputs%eqfile)) open(newunit = funit, file = self%filename, status='old',form='formatted') read(funit, '(5(a10),i2,i4,i4)') char, char, char, char, char, i, nw_in, nh_in allocate(spsi_bar(nw_in), sefit_R(nw_in), sefit_Z(nh_in), dummy(nw_in)) allocate(sefit_psi(nw_in, nh_in), f(nw_in), p(nw_in), q(nw_in)) read(funit, 2020) rwid, zhei, rcentr, rleft, zmid read(funit, 2020) self%R_mag, self%Z_mag, self%psi_0, self%psi_a, bcentr read(funit, 2020) xdum,xdum,xdum,xdum,xdum read(funit, 2020) xdum,xdum,xdum,xdum,xdum read(funit, 2020) (f(j), j = 1, nw_in) read(funit, 2020) (p(j), j = 1, nw_in) read(funit, 2020) (dummy(j), j = 1, nw_in) read(funit, 2020) (dummy(j), j = 1 ,nw_in) read(funit, 2020) ((sefit_psi(i,j) , i = 1, nw_in) , j = 1, nh_in) read(funit, 2020) (q(j), j = 1, nw_in) read(funit, '(2i5)') nbbbs, ndum allocate(rbbbs(nbbbs), zbbbs(nbbbs)) read(funit, 2020) (rbbbs(i), zbbbs(i) , i = 1, nbbbs) close(funit) ! pbar is defined by pbar = (psi-psi_0)/(psi_a-psi_0) fp and q are ! functions of pbar do i = 1, nw_in spsi_bar(i) = float(i - 1) / float(nw_in - 1) sefit_R(i) = rleft + rwid * float(i - 1) / float(nw_in - 1) end do do j = 1, nh_in sefit_Z(j) = ((float(j - 1) / float(nh_in - 1) ) - 0.5) * zhei !Should we add zmid? end do call self%shared_setup(nw_in, nh_in, spsi_bar, sefit_R, sefit_Z, f, p, q, sefit_psi, & nbbbs, rbbbs, zbbbs, bcentr, rleft, rwid, -0.5*zhei, zhei, inputs%big, .true.) !Should -0.5*zhei be zmid 2020 format (5e16.9) end subroutine efit_read_and_set !> FIXME : Add documentation pure subroutine dealloc_special_arrays(self) implicit none class(eeq_base_type), intent(in out) :: self deallocate(self%thetab, self%r_bound, self%eq_R, self%eq_Z) end subroutine dealloc_special_arrays subroutine setup_special_splines(self) use splines, only: new_periodic_spline use constants, only: twopi implicit none class(eeq_base_type), intent(in out) :: self self%rbound_spl = new_periodic_spline(self%thetab, self%r_bound, twopi, & drop_last_point = .true.) end subroutine setup_special_splines subroutine delete_special_splines(self) use splines, only: delete_periodic_spline implicit none class(eeq_base_type), intent(in out) :: self call delete_periodic_spline(self%rbound_spl) end subroutine delete_special_splines !> FIXME : Add documentation real function rcenter (self, rp) use constants, only: pi implicit none class(eeq_base_type), intent(in) :: self real, intent (in) :: rp rcenter = self%R_mag + 0.5 * (self%rfun(rp, 0.) - self%rfun(rp, pi)) end function rcenter end module eeq