shared_setup_cart Subroutine

public subroutine shared_setup_cart(self, nw_in, nh_in, spsi_bar, sefit_R, sefit_Z, f, p, q, sefit_psi, nbbbs, rbbbs, zbbbs, bcentr, rleft, rwid, zoffset, zhei, big, calc_aminor)

Type Bound

abstract_eqfile_cart_geo_type

Arguments

Type IntentOptional Attributes Name
class(abstract_eqfile_cart_geo_type), intent(inout) :: self
integer, intent(in) :: nw_in
integer, intent(in) :: nh_in
real, intent(in), dimension(:) :: spsi_bar
real, intent(in), dimension(:) :: sefit_R
real, intent(in), dimension(:) :: sefit_Z
real, intent(in), dimension(:) :: f
real, intent(in), dimension(:) :: p
real, intent(in), dimension(:) :: q
real, intent(in), dimension(:, :) :: sefit_psi
integer, intent(in) :: nbbbs
real, intent(inout), dimension(:) :: rbbbs
real, intent(inout), dimension(:) :: zbbbs
real, intent(in) :: bcentr
real, intent(in) :: rleft
real, intent(in) :: rwid
real, intent(in) :: zoffset
real, intent(in) :: zhei
integer, intent(in) :: big
logical, intent(in) :: calc_aminor

Contents

Source Code


Source Code

  subroutine shared_setup_cart(self, nw_in, nh_in, spsi_bar, sefit_R, sefit_Z, f, p, q, &
       sefit_psi, nbbbs, rbbbs, zbbbs, bcentr, rleft, rwid, zoffset, zhei, big, calc_aminor)
    use splines, only: inter_cspl, fitp_surf2, fitp_surf1
    use constants, only: twopi
    implicit none
    class(abstract_eqfile_cart_geo_type), intent(in out) :: self
    integer, intent(in) :: nw_in, nh_in, nbbbs, big
    real, dimension(:), intent(in) :: spsi_bar, sefit_R, sefit_Z, f, p, q
    real, dimension(:), intent(in out) :: rbbbs, zbbbs
    real, intent(in) :: bcentr, rleft, rwid, zoffset, zhei
    real, dimension(:, :), intent(in) :: sefit_psi
    logical, intent(in) :: calc_aminor
    real, dimension(:), allocatable :: zp, temp, zx1, zxm, zy1, zyn, eq_R, eq_Z
    real :: zxy11, zxym1, zxy1n, zxymn, psi_N, f_N
    integer :: i, j, ierr
    self%nw = nw_in * big ; self%nr = self%nw
    self%nh = nh_in * big ; self%nt = self%nh
    call self%alloc_arrays()
    do i = 1, self%nw
       self%psi_bar(i) = float(i - 1) / float(self%nw - 1)
    end do
    call inter_cspl(spsi_bar, f, self%psi_bar, self%fp)
    call inter_cspl(spsi_bar, p, self%psi_bar, self%pressure)
    call inter_cspl(spsi_bar, q, self%psi_bar, self%qsf)

    allocate(eq_R(self%nw), eq_Z(self%nh))
    do i = 1, self%nw
       eq_R(i) = rleft + rwid * float(i - 1) / float(self%nw - 1)
    end do

    do j = 1, self%nh
       eq_Z(j) = zoffset + ((float(j - 1) / float(self%nh - 1) )) * zhei
    end do

    ! Setup 2D interpolation of psi(R, Z)
    allocate(zp(3 * nw_in * nh_in), temp(nw_in + 2 * nh_in))
    allocate(zx1(nh_in), zxm(nh_in), zy1(nw_in), zyn(nw_in))
    call fitp_surf1(nw_in, nh_in, sefit_R, sefit_Z, sefit_psi, &
         nw_in, zx1, zxm, zy1, zyn, zxy11, zxym1, zxy1n, zxymn, &
         255, zp, temp, 1., ierr)

    do j = 1, self%nh
       do i = 1, self%nw
          self%eqpsi_2d(i, j) = fitp_surf2(eq_R(i), eq_Z(j), nw_in, nh_in, &
               sefit_R, sefit_Z, sefit_psi, nw_in, zp, 1.)
       end do
    end do
    deallocate(zp, temp)

    allocate(self%thetab(nbbbs), self%r_bound(nbbbs))
    self%thetab = atan2((zbbbs - self%Z_mag), (rbbbs - self%R_mag))
    self%r_bound = hypot((rbbbs - self%R_mag), (zbbbs - self%Z_mag))
    call sort(self%thetab, self%r_bound, zbbbs, rbbbs)

    ! Allow for duplicated points near +- pi:
    if (self%thetab(1) == self%thetab(2)) then
       self%thetab(1) = self%thetab(1) + twopi
       call sort(self%thetab, self%r_bound, zbbbs, rbbbs)
    end if

    if (self%thetab(nbbbs-1) == self%thetab(nbbbs)) then
       self%thetab(nbbbs) = self%thetab(nbbbs) - twopi
       call sort(self%thetab, self%r_bound, zbbbs, rbbbs)
    end if

    ! It isn't likely that a duplicate point would exist near theta = 0,
    ! so I am not allowing this possibility for now.
    do i = 1, nbbbs-1
       if (self%thetab(i+1) == self%thetab(i)) then
          ! put in kluge for duplicate points, which happens near theta=0:
          self%thetab(i+1) = self%thetab(i+1) +  1.e-8
       end if
    end do

    if (calc_aminor) call a_minor(rbbbs, zbbbs, self%Z_mag, self%aminor)

    self%B_psi = 0.0 ; self%diam = 0.0 ; self%rc = 0.0
    self%r_bound = self%r_bound / self%aminor
    self%R_mag = self%R_mag / self%aminor
    self%Z_mag = self%Z_mag / self%aminor
    eq_R = eq_R / self%aminor ; eq_Z = eq_Z / self%aminor
    self%eq_R = eq_R ; self%eq_Z = eq_Z
    do j = 1, self%nh
       self%R_psi(:, j) = eq_R ; self%Z_psi(:, j) = eq_Z(j)
    end do

    self%B_T = abs(bcentr)
    psi_N = self%B_T * self%aminor**2
    self%psi_a = self%psi_a / psi_N
    self%psi_0 = self%psi_0 / psi_N
    self%eqpsi_2d = self%eqpsi_2d / psi_N
    f_N = self%B_T * self%aminor
    self%fp = self%fp / f_N
    ! MKS: beta = 2 mu_0 p / B**2
    self%beta = 4. * twopi * self%pressure * 1.e-7 / self%B_T**2
    self%beta_0 = self%beta(1)
    self%pressure = self%pressure / self%pressure(1)

    do j = 1, self%nh
       do i = 1, self%nw
          if (eq_Z(j) == self%Z_mag .and. eq_R(i) == self%R_mag) then
             self%eqth(i, j) = 0.  ! value should not matter
          else
             self%eqth(i, j) = atan2( (eq_Z(j) - self%Z_mag), (eq_R(i) - self%R_mag))
          end if
       end do
    end do
    self%eqpsi = self%psi_bar * (self%psi_a - self%psi_0) + self%psi_0
    self%has_full_theta_range = .true.
    if (.true.) then
       write (*,*) "Finished shared_setup... imported ",self%type_name," equilibrium"
       write (*,*) 'Some important quantities:'
       write (*,*) "aminor", self%aminor
       write (*,*) 'R_mag', self%R_mag
       write (*,*) 'B_T0', self%B_T
       write (*,*) 'beta', self%beta_0
    end if
  end subroutine shared_setup_cart