module slab_equilibrium_m !! Definition of slab equilibrium !! N.b. this equilibrium is axisymmetric - phi, when required, is not used use precision_m, only: FP, FP_EPS use screen_io_m, only: get_stdout use descriptors_m, only: DISTRICT_WALL, DISTRICT_SOL, & DISTRICT_CLOSED, DISTRICT_CORE use equilibrium_m, only: equilibrium_t use params_equi_slab_m, only: read_params_slab, get_slab_boxsize, & get_slab_sol, get_slab_yperiodic implicit none type, public, extends(equilibrium_t) :: slab_t !! Type that represents the slab equilibrium logical, private :: sol !! Indicates if slab is open or closed field lines logical, private :: yperiodic !! Indicates if periodic in y direction contains procedure, public :: init procedure, public :: display procedure, public :: debug procedure, public :: is_axisymmetric procedure, public :: rho procedure, public :: bx procedure, public :: by procedure, public :: btor procedure, public :: jacobian procedure, public :: epol procedure, public :: erad procedure, public :: district procedure, public :: in_vessel procedure, public :: mag_axis_loc procedure, public :: is_yperiodic end type contains subroutine init(self, filename) !! Initialises the slab equilibrium class(slab_t), intent(inout) :: self !! Instance of the type character(len=*), intent(in), optional :: filename !! Filename, where parameters are read from if (present(filename)) then call read_params_slab(filename) endif self%x0 = 0.0_FP self%y0 = 0.0_FP self%xmin = 0.0_FP self%xmax = get_slab_boxsize() self%ymin = 0.0_FP self%ymax = get_slab_boxsize() self%rhomin = self%xmin self%rhomax = self%xmax self%sol = get_slab_sol() self%yperiodic = get_slab_yperiodic() self%initialized = .true. end subroutine init subroutine display(self) !! Prints to console information about the equilibrium class(slab_t), intent(in) :: self write(get_stdout(),99) & self%sol, & self%yperiodic 99 FORMAT( / & 1X,'Slab equilibrium:'/, & 1X,'sol = ',L1 /, & 1X,'yperiodic = ',L1 & ) end subroutine display subroutine debug(self) !! Prints to console information about the equilibrium class(slab_t), intent(in) :: self write(get_stdout(),100) & self%x0, & self%y0, & self%rhomin, & self%rhomax, & self%xmin, & self%xmax, & self%ymin, & self%ymax, & self%sol, & self%yperiodic 100 FORMAT( / & 1X,'Slab equilibrium:'/, & 1X,'x0 = ',ES14.6E3 /, & 1X,'y0 = ',ES14.6E3 /, & 1X,'rhomin = ',ES14.6E3 /, & 1X,'rhomax = ',ES14.6E3 /, & 1X,'xmin = ',ES14.6E3 /, & 1X,'xmax = ',ES14.6E3 /, & 1X,'ymin = ',ES14.6E3 /, & 1X,'ymax = ',ES14.6E3 /, & 1X,'sol = ',L1 /, & 1X,'yperiodic = ',L1 & ) end subroutine debug logical function is_axisymmetric(self) class(slab_t), intent(in) :: self is_axisymmetric = .true. end function real(FP) function rho(self, x, y, phi) !! flux surface label = x class(slab_t),intent(in) :: self real(FP), intent(in) :: x, y, phi rho = x end function rho real(FP) function bx(self, x, y, phi) !! magnetic field component bx = 0 class(slab_t),intent(in) :: self real(FP), intent(in) :: x, y, phi bx = 0.0_FP end function bx real(FP) function by(self, x, y, phi) !! magnetic field component by = 0 class(slab_t),intent(in) :: self real(FP), intent(in) :: x, y, phi by = 0.0_FP end function by real(FP) function btor(self, x, y, phi) !! magnetic field strength normalised = 1 class(slab_t),intent(in) :: self real(FP), intent(in) :: x, y, phi btor = 1.0_FP end function btor real(FP) function jacobian(self, x, y, phi) !! Jacobian of geometry (J=1) class(slab_t),intent(in) :: self real(FP), intent(in) :: x, y, phi jacobian = 1.0_FP end function jacobian subroutine epol(self, x, y, phi, epolx, epoly) !! unit vector along poloidal (along flux surface) direction class(slab_t),intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP), intent(out) :: epolx !! unit vector component in x-direction real(FP), intent(out) :: epoly !! unit vector component in y-direction epolx = 0.0_FP epoly = 1.0_FP end subroutine epol subroutine erad(self, x, y, phi, eradx, erady) !! unit vector along radial (across flux surface) direction class(slab_t),intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP), intent(out) :: eradx !! unit vector component in x-direction real(FP), intent(out) :: erady !! unit vector component in y-direction eradx = 1.0_FP erady = 0.0_FP end subroutine erad integer function district(self, x, y, phi) !! returns in which district point (x, y, phi) is (see module descriptors_m) class(slab_t), intent(in) :: self real(FP), intent(in) :: x, y, phi if (self%sol) then if (x <= self%xmin + FP_EPS) then district = DISTRICT_WALL elseif (x >= self%xmax - FP_EPS) then district = DISTRICT_WALL else district = DISTRICT_SOL endif else if (x <= self%xmin + FP_EPS) then district = DISTRICT_CORE elseif (x >= self%xmax - FP_EPS) then district = DISTRICT_WALL else district = DISTRICT_CLOSED endif endif if (.not. self%yperiodic) then if (y <= self%ymin + FP_EPS) then district = DISTRICT_WALL elseif (y >= self%ymax - FP_EPS) then district = DISTRICT_WALL endif endif end function district logical function in_vessel(self, x, y, phi) !! For slab always inside vessel class(slab_t), intent(in) :: self real(FP), intent(in) :: x, y, phi in_vessel = .true. end function in_vessel subroutine mag_axis_loc(self, phi, axis_x, axis_y) !! Returns the coordinates of magnetic axis class(slab_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: phi !! Toroidal angle real(FP), intent(out) :: axis_x !! x-coordinate of the magnetic axis real(FP), intent(out) :: axis_y !! y-coordinate of the magnetic axis axis_x = self%x0 axis_y = self%y0 end subroutine mag_axis_loc pure logical function is_yperiodic(self) !! Returns .true. if the equilibrium is y periodic and .false. !! otherwise class(slab_t), intent(in) :: self is_yperiodic = self%yperiodic end function is_yperiodic end module slab_equilibrium_m