read_and_set Subroutine

private subroutine read_and_set(self, inputs)

Uses

This subroutine reads an DFIT output file containing the axisymmetric magnetic field geometry on a rectangular domain defined by the coordinates (R,Z).

Type Bound

deq_type

Arguments

Type IntentOptional Attributes Name
class(deq_type), intent(inout) :: self
type(geo_input_type), intent(in) :: inputs

Contents

Source Code


Source Code

  subroutine read_and_set(self, inputs)
    use mp, only: mp_abort
    use geo_utils, only: sort
    implicit none
    class(deq_type), intent(in out) :: self
    type(geo_input_type), intent(in) :: inputs
    real, dimension(:), allocatable :: spsi_bar, sdfit_R, sdfit_Z, rbbbs, zbbbs, eq_R, eq_Z, f, p, q
    real, dimension(:, :), allocatable :: sdfit_psi
    real :: rwid, zhei, bcentr, delta_R, delta_Z, lref, bref
    integer :: i, j, jmin, jmax, nw_in, nh_in, funit, nbbbs, nrho
    logical :: eqnorm_exists
    self%type_name = 'deq'
    ! Need to generalize initialization condition if equilibrium changes
    lref = 0. ; bref = 0.
    ! read eqnorm file if exists
    inquire(file = trim(inputs%eqnormfile), exist = eqnorm_exists)
    if (eqnorm_exists) then
       open(newunit = funit, file = trim(inputs%eqnormfile), status='old', form='formatted')
       read(funit, *) lref, bref
       close(funit)
    end if

    self%filename = trim(adjustl(inputs%eqfile))
    open(newunit = funit, file = self%filename, status='old', form='formatted')

    ! Read the data
    read(funit, *) nw_in, nh_in, nrho
    allocate (self%rho_mid(nrho), self%psi_mid(nrho))
    allocate(spsi_bar(nw_in), sdfit_R(nw_in), sdfit_Z(nh_in), sdfit_psi(nw_in, nh_in))
    allocate(eq_R(self%nw), eq_Z(self%nh), f(nw_in), p(nw_in), q(nw_in))
    read(funit, 2020) rwid, zhei
    read(funit, 2020) self%R_mag, self%Z_mag
    read(funit, 2020) self%psi_0, self%psi_a, bcentr
    f = 0 ; q = 0
    read(funit, 2020) (p(j), j = 1, nw_in)  ! pressure read
    read(funit, 2020) ((sdfit_psi(i,j) , i = 1, nw_in) , j = 1, nh_in)

    ! Only needed for irho = 4
    read (funit, 2020) (self%rho_mid(i), i = 1, nrho)
    read (funit, 2020) (self%psi_mid(i), i = 1, nrho)
    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)
       sdfit_R(i) = rwid * float(i - 1) / float(nw_in - 1)
    end do

    do j = 1, nh_in
       sdfit_Z(j) = (float(j - 1) / float(nh_in - 1) - 0.5) * zhei
    end do

    ! stay away from edge of box
    delta_R = 1.02 ; delta_Z = 1.02

    nbbbs = 2 * (nw_in + nh_in) - 5
    allocate(rbbbs(nbbbs), zbbbs(nbbbs))
    jmin = 1 ; jmax = nh_in / 2
    do j = jmin, jmax
       rbbbs(j) = (sdfit_R(1) - 0.5 * rwid) / delta_R + 0.5 * rwid
       zbbbs(j) = sdfit_Z(j + nh_in / 2) / delta_Z
    end do

    jmin = jmax + 1 ; jmax = jmax + nw_in - 2
    do j = jmin, jmax
       rbbbs(j) = (sdfit_R(j - jmin + 2) - 0.5 * rwid) / delta_R + 0.5 * rwid
       zbbbs(j) = sdfit_Z(nh_in) / delta_Z
    end do

    jmin = jmax + 1 ; jmax = jmax + nh_in - 1
    do j = jmin, jmax
       rbbbs(j) = (sdfit_R(nw_in) - 0.5 * rwid) / delta_R + 0.5 * rwid
       zbbbs(j) = sdfit_Z(jmax - j + 2) / delta_Z
    end do

    jmin = jmax + 1 ; jmax = jmax + nw_in - 1
    do j = jmin, jmax
       rbbbs(j) = (sdfit_R(jmax - j + 2) - 0.5 * rwid) / delta_R + 0.5 * rwid
       zbbbs(j) = sdfit_Z(1) / delta_Z
    end do

    jmin = jmax + 1 ; jmax = jmax + nh_in / 2
    do j = jmin, jmax
       rbbbs(j) = (sdfit_R(1) - 0.5 * rwid) / delta_R + 0.5 * rwid
       zbbbs(j) = sdfit_Z(j - jmin + 2) / delta_Z
    end do

    do j = 1, nbbbs-1
       if ((rbbbs(j) == rbbbs(j+1)) .and. (zbbbs(j) == zbbbs(j+1))) then
          write(*,*) 'duplicates: ',rbbbs(j),zbbbs(j)
       end if
    end do

    ! Overwrite inputs with normalisations if provided
    if (bref /= 0.) then
       bcentr = bref
       write(6,*) 'B_T0 = ', bref
    end if

    if (lref /= 0.) then
       self%aminor = lref
       write(6,*) 'aminor = ', self%aminor, lref
    else
       self%aminor = self%R_mag
    end if

    call self%shared_setup(nw_in, nh_in, spsi_bar, sdfit_R, sdfit_Z, f, p, q, sdfit_psi, &
         nbbbs, rbbbs, zbbbs, bcentr, 0.0, rwid, -0.5*zhei, zhei, 1, .false.)
2020 format (5e16.8)
  end subroutine read_and_set