mod_equilibrium.f90

Go to the documentation of this file.

module mod_equilibrium

  use itm_types
  use mod_parameter
  use euitm_schemas
  use mod_helena_io

  implicit none

  public :: equilibrium_constructor, equilibrium_destructor
  public :: fill_equilibrium

  contains

!--------------------------------------------------------------------------
  subroutine equilibrium_constructor(equilibrium, nr, nchi)
!--------------------------------------------------------------------------

    implicit none

    integer(itm_i4) :: nr, nchi
    type (type_equilibrium) :: equilibrium

!-- allocate pointers

!-- 1d profiles
    allocate(equilibrium%profiles_1d%psi(nr))
    allocate(equilibrium%profiles_1d%phi(nr))
    allocate(equilibrium%profiles_1d%pressure(nr))
    allocate(equilibrium%profiles_1d%F_dia(nr))
    allocate(equilibrium%profiles_1d%pprime(nr))
    allocate(equilibrium%profiles_1d%ffprime(nr))
    allocate(equilibrium%profiles_1d%jphi(nr))
    allocate(equilibrium%profiles_1d%jparallel(nr))
    allocate(equilibrium%profiles_1d%q(nr))
    allocate(equilibrium%profiles_1d%r_inboard(nr))
    allocate(equilibrium%profiles_1d%r_outboard(nr))
    allocate(equilibrium%profiles_1d%rho_tor(nr))
    allocate(equilibrium%profiles_1d%dpsidrho_tor(nr))
    allocate(equilibrium%profiles_1d%rho_vol(nr))
    allocate(equilibrium%profiles_1d%beta_pol(nr))
    allocate(equilibrium%profiles_1d%li(nr))
    allocate(equilibrium%profiles_1d%elongation(nr))
    allocate(equilibrium%profiles_1d%tria_upper(nr))
    allocate(equilibrium%profiles_1d%tria_lower(nr))
    allocate(equilibrium%profiles_1d%volume(nr))
    allocate(equilibrium%profiles_1d%vprime(nr))
    allocate(equilibrium%profiles_1d%dvdrho(nr))
    allocate(equilibrium%profiles_1d%area(nr))
    allocate(equilibrium%profiles_1d%aprime(nr))
    allocate(equilibrium%profiles_1d%surface(nr))
    allocate(equilibrium%profiles_1d%ftrap(nr))
    allocate(equilibrium%profiles_1d%gm1(nr))
    allocate(equilibrium%profiles_1d%gm2(nr))
    allocate(equilibrium%profiles_1d%gm3(nr))
    allocate(equilibrium%profiles_1d%gm4(nr))
    allocate(equilibrium%profiles_1d%gm5(nr))
    allocate(equilibrium%profiles_1d%gm6(nr))
    allocate(equilibrium%profiles_1d%gm7(nr))
    allocate(equilibrium%profiles_1d%gm8(nr))
    allocate(equilibrium%profiles_1d%gm9(nr))
    allocate(equilibrium%profiles_1d%b_av(nr))
    allocate(equilibrium%profiles_1d%b_min(nr))
    allocate(equilibrium%profiles_1d%b_max(nr))

!-- 2d profiles
    allocate(equilibrium%profiles_2d(1))
    allocate(equilibrium%profiles_2d(1)%grid_type(4))
    allocate(equilibrium%profiles_2d(1)%grid%dim1(nr))
    allocate(equilibrium%profiles_2d(1)%grid%dim2(nchi))
    allocate(equilibrium%profiles_2d(1)%r(nr, nchi))
    allocate(equilibrium%profiles_2d(1)%z(nr, nchi))
    allocate(equilibrium%profiles_2d(1)%psi(nr, nchi))
    allocate(equilibrium%profiles_2d(1)%theta(nr, nchi))
    allocate(equilibrium%profiles_2d(1)%jphi(nr, nchi))
    allocate(equilibrium%profiles_2d(1)%jpar(nr, nchi))
    allocate(equilibrium%profiles_2d(1)%br(nr, nchi))
    allocate(equilibrium%profiles_2d(1)%bz(nr, nchi))
    allocate(equilibrium%profiles_2d(1)%bphi(nr, nchi))

!-- boundary
    allocate(equilibrium%eqgeometry%boundary(1))
    allocate(equilibrium%eqgeometry%boundary(1)%r(nchi))
    allocate(equilibrium%eqgeometry%boundary(1)%z(nchi))

!-- grid
    allocate(equilibrium%coord_sys%grid_type(4))
    allocate(equilibrium%coord_sys%grid%dim1(nr))
    allocate(equilibrium%coord_sys%grid%dim2(nchi))
    allocate(equilibrium%coord_sys%position%r(nr, nchi))
    allocate(equilibrium%coord_sys%position%z(nr, nchi))
    allocate(equilibrium%coord_sys%g_11(nr, nchi))
    allocate(equilibrium%coord_sys%g_12(nr, nchi))
    allocate(equilibrium%coord_sys%g_13(nr, nchi))
    allocate(equilibrium%coord_sys%g_22(nr, nchi))
    allocate(equilibrium%coord_sys%g_23(nr, nchi))
    allocate(equilibrium%coord_sys%g_33(nr, nchi))
    allocate(equilibrium%coord_sys%jacobian(nr, nchi))

!-- initialize pointers
    equilibrium%profiles_1d%psi = 0._r8
    equilibrium%profiles_1d%phi = 0._r8
    equilibrium%profiles_1d%pressure = 0._r8
    equilibrium%profiles_1d%F_dia = 0._r8
    equilibrium%profiles_1d%pprime = 0._r8
    equilibrium%profiles_1d%ffprime = 0._r8
    equilibrium%profiles_1d%jphi = 0._r8
    equilibrium%profiles_1d%jparallel = 0._r8
    equilibrium%profiles_1d%q = 0._r8
    equilibrium%profiles_1d%r_inboard = 0._r8
    equilibrium%profiles_1d%r_outboard = 0._r8
    equilibrium%profiles_1d%rho_tor = 0._r8
    equilibrium%profiles_1d%dpsidrho_tor = 0._r8
    equilibrium%profiles_1d%rho_vol = 0._r8
    equilibrium%profiles_1d%beta_pol = 0._r8
    equilibrium%profiles_1d%li = 0._r8
    equilibrium%profiles_1d%elongation = 0._r8
    equilibrium%profiles_1d%tria_upper = 0._r8
    equilibrium%profiles_1d%tria_lower = 0._r8
    equilibrium%profiles_1d%volume = 0._r8
    equilibrium%profiles_1d%vprime = 0._r8
    equilibrium%profiles_1d%dvdrho = 0._r8
    equilibrium%profiles_1d%area = 0._r8
    equilibrium%profiles_1d%aprime = 0._r8
    equilibrium%profiles_1d%surface = 0._r8
    equilibrium%profiles_1d%ftrap = 0._r8
    equilibrium%profiles_1d%gm1 = 0._r8
    equilibrium%profiles_1d%gm2 = 0._r8
    equilibrium%profiles_1d%gm3 = 0._r8
    equilibrium%profiles_1d%gm4 = 0._r8
    equilibrium%profiles_1d%gm5 = 0._r8
    equilibrium%profiles_1d%gm6 = 0._r8
    equilibrium%profiles_1d%gm7 = 0._r8
    equilibrium%profiles_1d%gm8 = 0._r8
    equilibrium%profiles_1d%gm9 = 0._r8
    equilibrium%profiles_1d%b_av = 0._r8
    equilibrium%profiles_1d%b_min = 0._r8
    equilibrium%profiles_1d%b_max = 0._r8

    equilibrium%profiles_2d(1)%grid%dim1 = 0._r8
    equilibrium%profiles_2d(1)%grid%dim2 = 0._r8
    equilibrium%profiles_2d(1)%r = 0._r8
    equilibrium%profiles_2d(1)%z = 0._r8
    equilibrium%profiles_2d(1)%psi = 0._r8
    equilibrium%profiles_2d(1)%theta = 0._r8
    equilibrium%profiles_2d(1)%jphi = 0._r8
    equilibrium%profiles_2d(1)%jpar = 0._r8
    equilibrium%profiles_2d(1)%br = 0._r8
    equilibrium%profiles_2d(1)%bz = 0._r8
    equilibrium%profiles_2d(1)%bphi = 0._r8

    equilibrium%eqgeometry%boundary(1)%r = 0._r8
    equilibrium%eqgeometry%boundary(1)%z = 0._r8

    equilibrium%coord_sys%grid%dim1 = 0._r8
    equilibrium%coord_sys%grid%dim2 = 0._r8
    equilibrium%coord_sys%position%r = 0._r8
    equilibrium%coord_sys%position%z = 0._r8
    equilibrium%coord_sys%g_11 = 0._r8
    equilibrium%coord_sys%g_12 = 0._r8
    equilibrium%coord_sys%g_13 = 0._r8
    equilibrium%coord_sys%g_22 = 0._r8
    equilibrium%coord_sys%g_23 = 0._r8
    equilibrium%coord_sys%g_33 = 0._r8
    equilibrium%coord_sys%jacobian = 0._r8

  end subroutine equilibrium_constructor

!--------------------------------------------------------------------------
  subroutine equilibrium_destructor(equilibrium)
!--------------------------------------------------------------------------

    implicit none

    type (type_equilibrium) :: equilibrium

!-- deallocate pointers

!-- 1d profiles
    deallocate(equilibrium%profiles_1d%psi)
    deallocate(equilibrium%profiles_1d%phi)
    deallocate(equilibrium%profiles_1d%pressure)
    deallocate(equilibrium%profiles_1d%F_dia)
    deallocate(equilibrium%profiles_1d%pprime)
    deallocate(equilibrium%profiles_1d%ffprime)
    deallocate(equilibrium%profiles_1d%jphi)
    deallocate(equilibrium%profiles_1d%jparallel)
    deallocate(equilibrium%profiles_1d%q)
    deallocate(equilibrium%profiles_1d%r_inboard)
    deallocate(equilibrium%profiles_1d%r_outboard)
    deallocate(equilibrium%profiles_1d%rho_tor)
    deallocate(equilibrium%profiles_1d%dpsidrho_tor)
    deallocate(equilibrium%profiles_1d%rho_vol)
    deallocate(equilibrium%profiles_1d%beta_pol)
    deallocate(equilibrium%profiles_1d%li)
    deallocate(equilibrium%profiles_1d%elongation)
    deallocate(equilibrium%profiles_1d%tria_upper)
    deallocate(equilibrium%profiles_1d%tria_lower)
    deallocate(equilibrium%profiles_1d%volume)
    deallocate(equilibrium%profiles_1d%vprime)
    deallocate(equilibrium%profiles_1d%dvdrho)
    deallocate(equilibrium%profiles_1d%area)
    deallocate(equilibrium%profiles_1d%aprime)
    deallocate(equilibrium%profiles_1d%surface)
    deallocate(equilibrium%profiles_1d%ftrap)
    deallocate(equilibrium%profiles_1d%gm1)
    deallocate(equilibrium%profiles_1d%gm2)
    deallocate(equilibrium%profiles_1d%gm3)
    deallocate(equilibrium%profiles_1d%gm4)
    deallocate(equilibrium%profiles_1d%gm5)
    deallocate(equilibrium%profiles_1d%gm6)
    deallocate(equilibrium%profiles_1d%gm7)
    deallocate(equilibrium%profiles_1d%gm8)
    deallocate(equilibrium%profiles_1d%gm9)
    deallocate(equilibrium%profiles_1d%b_av)
    deallocate(equilibrium%profiles_1d%b_min)
    deallocate(equilibrium%profiles_1d%b_max)

!-- 2d profiles
    deallocate(equilibrium%profiles_2d(1)%grid_type)
    deallocate(equilibrium%profiles_2d(1)%grid%dim1)
    deallocate(equilibrium%profiles_2d(1)%grid%dim2)
    deallocate(equilibrium%profiles_2d(1)%r)
    deallocate(equilibrium%profiles_2d(1)%z)
    deallocate(equilibrium%profiles_2d(1)%psi)
    deallocate(equilibrium%profiles_2d(1)%theta)
    deallocate(equilibrium%profiles_2d(1)%jphi)
    deallocate(equilibrium%profiles_2d(1)%jpar)
    deallocate(equilibrium%profiles_2d(1)%br)
    deallocate(equilibrium%profiles_2d(1)%bz)
    deallocate(equilibrium%profiles_2d(1)%bphi)
    deallocate(equilibrium%profiles_2d)

!-- boundary
    deallocate(equilibrium%eqgeometry%boundary(1)%r)
    deallocate(equilibrium%eqgeometry%boundary(1)%z)
    deallocate(equilibrium%eqgeometry%boundary)

!-- grid
    deallocate(equilibrium%coord_sys%grid_type)
    deallocate(equilibrium%coord_sys%grid%dim1)
    deallocate(equilibrium%coord_sys%grid%dim2)
    deallocate(equilibrium%coord_sys%position%r)
    deallocate(equilibrium%coord_sys%position%z)
    deallocate(equilibrium%coord_sys%g_11)
    deallocate(equilibrium%coord_sys%g_12)
    deallocate(equilibrium%coord_sys%g_13)
    deallocate(equilibrium%coord_sys%g_22)
    deallocate(equilibrium%coord_sys%g_23)
    deallocate(equilibrium%coord_sys%g_33)
    deallocate(equilibrium%coord_sys%jacobian)

  end subroutine equilibrium_destructor

!-----------------------------------------------------------------------------
  subroutine fill_equilibrium(a, xaxis, yaxis, current_averaging, equilibrium)
!-----------------------------------------------------------------------------
! This subroutine fills the equilibrium CPO equilibrium through calls
! of the appropriate subroutines.
! This subroutine requires equilibrium to be constructed by a
! call to equilibrium_constructor and the mesh to be aligned with
! the flux surfaces (possibly call to remesh necessary).
!-----------------------------------------------------------------------------
    use mod_dat, only : eps, rvac, bvac, alfa, hbt, b, verbosity, p_sep, &
                        bmag, pscale, r_vessel, zgeo, cpsurfin
    use mod_nodes, only : psikn
    use mod_profiles
    use mod_spline
    use current_calculations
    use global_plasma_parameters

    implicit none

    type (type_equilibrium), intent(inout) :: equilibrium
    character(len = *), intent(in) :: current_averaging
    real(r8), intent(in) :: a, xaxis, yaxis

    real(r8) :: p_out(nr), rbphi_out(nr), dp_out(nr), drbphi_out(nr)
    real(r8) :: bm
    real(r8) :: r0, b0
    real(r8) :: rbscale
    real(r8) :: raxis, rminor, rmag
    real(r8) :: cpsurf_out
    real(r8) :: toroidal_current
    integer(itm_i4) :: i

    real(r8) :: dp_dpsi, dgamma_dpsi

!-- set toroid_field, required for calculation of rho_tor
    equilibrium%global_param%toroid_field%r0 = r_vessel
    equilibrium%global_param%toroid_field%b0 = rvac * bvac / r_vessel

!-- profiles yields the pressure normalized to eps*Bvac^2/(mu_0 * alfa^2)
!   and R Bphi normalized to Rvac*Bvac (dp_out and drbphi_out are only
!   dummies here, since profiles gives their derivatives with respect to rho
    call profiles(p_out, rbphi_out, dp_out, drbphi_out, a)
!-- caculate dp_hat/dpsi_hat and drbphi_hat/dpsi_hat
    do i = 1, nr
      if (hbt) then
        dp_out(i) = -0.5_r8 * a * b * dp_dpsi(psikn(nr - i + 1))
        drbphi_out(i) = dp_out(i) - eps * a * dgamma_dpsi(psikn(nr - i + 1))
      else
        dp_out(i) = -eps * a * b * dp_dpsi(psikn(nr - i + 1))
        drbphi_out(i) = -eps * a * dgamma_dpsi(psikn(nr - i + 1))
      end if
      drbphi_out(i) = 1._r8 / (2._r8 * rbphi_out(i)) * 2._r8 * eps &
       * drbphi_out(i) / alfa**2
    end do
    raxis = 1._r8 + eps * xaxis
!-- normalized Bmag/Bvac on most recent axis
    bm = rbphi_out(1) / raxis
!-- Bmag on most recent axis
    bmag = bm * bvac
!-- change normalization to Bmag**2/mu_0 for p_out and current Rmag * Bmag
!   derivatives with respect to psi_hat
    b0 = 1._r8/ bm
    r0 = 1._r8/ raxis
    pscale = b0**2 * eps / alfa**2
    rbscale = b0 * r0
    p_out = p_out * pscale
    dp_out = dp_out * pscale
    rbphi_out = rbphi_out * rbscale
    drbphi_out = drbphi_out * rbscale

!-- unnormalized poloidal flux psi [Vs] calculated from alfa
    cpsurf_out = twopi * eps**2 * rvac**2 * bvac / alfa  

    if (verbosity > 2) write(iu6, *) 'cpsurf_out = ', cpsurf_out

!-- in SI units
    rminor = eps * rvac
    rmag = raxis * rvac

!-- denormalize p_out, dp_out, rbphi_out, and drbphi_out
    p_out = p_out * bmag**2 / mu0
    dp_out = dp_out * bmag**2 / (mu0 * cpsurf_out)
    rbphi_out = rbphi_out * rmag * bmag
    drbphi_out = drbphi_out * rmag * bmag / cpsurf_out

!-- add pressure at separatrix since HELENA internally assumes p(sep) = 0
    p_out = p_out + p_sep

!-- calculate current densities and plasma current
    call current_calculations_constructor(nr, np)
    call total_current(a, xaxis, toroidal_current)
!-- re-calculate current_densities since total_current calculates current
!   densities with current_averaging = 'area'
    call current_densities(a, xaxis, current_averaging)
    call si_currents

!-- geometric and magnetic axis
    equilibrium%global_param%mag_axis%position%r = rmag
    equilibrium%global_param%mag_axis%position%z = zgeo + yaxis * rminor
    equilibrium%global_param%mag_axis%bphi = bmag
    equilibrium%eqgeometry%geom_axis%r = rvac
    equilibrium%eqgeometry%geom_axis%z = 0._r8   ! by definition
    equilibrium%eqgeometry%a_minor = rminor
!-- 1d profiles
    do i = 1, nr
      equilibrium%profiles_1d%psi(i) = cpsurf_out * psikn(nr + 1 - i)
    end do
    equilibrium%profiles_1d%psi(1) = 0._r8
    equilibrium%profiles_1d%pressure = p_out
    equilibrium%profiles_1d%pprime = dp_out
    equilibrium%profiles_1d%F_dia = rbphi_out
    equilibrium%profiles_1d%ffprime = drbphi_out * rbphi_out
    equilibrium%profiles_1d%jphi = j_tor_av
!-- cross section, surface area, and plasma volume
    equilibrium%global_param%area = cross_section(equilibrium, &
     rminor, eps)
    equilibrium%global_param%volume = plasma_volume(equilibrium, &
     rminor, eps)
    call surface_area(equilibrium, rminor, eps)
!-- plasma current
    equilibrium%global_param%i_plasma = i_tor(nr)
!-- calculate plasma beta and li
    call allocate_spline_coefficients(p_spline, npts)
    call spline(npts, psivec, p_int, dpres(1) * sign(1._r8, cpsurfin), &
     dpres(npts) * sign(1._r8, cpsurfin), 1, p_spline)
    equilibrium%global_param%beta_pol = beta_poloidal(a, equilibrium)
    equilibrium%global_param%beta_tor = beta_toroidal(a, equilibrium)
    equilibrium%global_param%li = internal_inductance(a, equilibrium)
!-- calculate W_MHD
    equilibrium%global_param%w_mhd = w_mhd(equilibrium, rminor, eps)
!-- Ip in MA, beta_normal in %
    equilibrium%global_param%beta_normal = 100._r8 * rminor &
     * abs(equilibrium%global_param%toroid_field%b0) &
     * equilibrium%global_param%beta_tor &
     / abs(equilibrium%global_param%i_plasma) * 1e6_r8
    call deallocate_spline_coefficients(p_spline)
!-- rho_vol following ITM definition
    equilibrium%profiles_1d%rho_vol = sqrt(equilibrium%profiles_1d%volume &
     / equilibrium%profiles_1d%volume(nr))

    if (verbosity > 2) &
     write(iu6, *) 'total current inside fill_equilibrium ', i_tor(nr)

    call current_calculations_destructor

  end subroutine fill_equilibrium

!-----------------------------------------------------------------------------
  subroutine adapt_alfa(a, xaxis, cx, cy, npts)
!-----------------------------------------------------------------------------
! This subroutine adapts the scaling parameter alfa, i.e., the poloidal
! flux cpsurf at the plasma boundary, using the plasma current Ip
! or q95 unless alfa itself is specified on input.
! A call to this subroutine is required after each change of the scaling
! parameter a which maintains the normalized poloidal flux at the boundary
! to be 1.
!-----------------------------------------------------------------------------

    use mod_dat, only : alfa, ip, q95, match, verbosity
    use mod_helena_io, only : iu6

    use current_calculations

    implicit none

    real(r8), intent(in) :: a, xaxis
    real(r8), intent(in) :: cx, cy
    integer(itm_i4), intent(in) :: npts

    real(r8) :: dummy(npts)
    real(r8) :: current
    real(r8) :: q95out, q1out

    if (trim(match) == "q95") then
      call safety_factor(xaxis, a, cx, cy, dummy, q95out, q1out)
      if (verbosity > 1) then
        write(iu6, *) 'recalculate alfa from safety factor q95 = ', q95
        write(iu6, *) 'q1out = ', q1out
        write(iu6, *) 'old alfa = ', alfa
      end if
      alfa = alfa * q95 / q1out
      if (verbosity > 1) then
        write(iu6, *) 'new alfa = ', alfa
      end if
      call total_current(a, xaxis, current)
      ip = current
    else if (trim(match) == "Ip") then
      call total_current(a, xaxis, current)
      if (verbosity > 1) then
        write(iu6, *) 'recalculate alfa from plasma current Ip = ', ip
        write(iu6, *) 'current Ip = ', ip
        write(iu6, *) 'new current = ', current
        write(iu6, *) 'old alfa = ', alfa
      end if
      alfa = current / ip
      if (verbosity > 1) write(iu6, *) 'new alfa = ', alfa
    else    ! match = "psi_bound"
      call total_current(a, xaxis, current)
    end if

  end subroutine adapt_alfa

end module mod_equilibrium