This subroutine reads a generic NetCDF equilibrium file containing the axisymmetric magnetic field geometry in flux coordinates
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(geq_type), | intent(inout) | :: | self | |||
| type(geo_input_type), | intent(in) | :: | inputs |
subroutine read_and_set(self, inputs)
implicit none
class(geq_type), intent(in out) :: self
type(geo_input_type), intent(in) :: inputs
# ifdef NETCDF
real, parameter :: psi_scale = 1.0e-8, b_scale = 1.0e-4, l_scale = 1.0e-2
real :: f_N, psi_N
integer :: istatus, ncid, id, j
self%type_name = 'geq'
self%filename = trim(adjustl(inputs%eqfile))
istatus = nf90_open(self%filename, NF90_NOWRITE, ncid)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, file = self%filename)
! netcdf read scalar: nr
!
! nr == number of radial grid points in radial eq grid
istatus = nf90_inq_dimid (ncid, 'psi', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, dim='psi')
istatus = nf90_inquire_dimension (ncid, id, len = self%nr)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, dimid=id)
! netcdf read scalar: nt
!
! nt == number of theta grid points in theta eq grid
istatus = nf90_inq_dimid (ncid, 'theta', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, dim='theta')
istatus = nf90_inquire_dimension (ncid, id, len = self%nt)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, dimid=id)
call self%alloc_arrays()
! netcdf read scalars: psi_0,psi_a,B_T,I_0
!
! psi_0 == value of psi at the magnetic axis
! psi_a == value of psi at the boundary of the plasma
! B_T == vacuum toroidal magnetic field at R_center
! I_0 == total plasma current
istatus = nf90_inq_varid (ncid, 'psi_0', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='psi_0')
istatus = nf90_get_var (ncid, id, self%psi_0)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
istatus = nf90_inq_varid (ncid, 'psi_a', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='psi_a')
istatus = nf90_get_var (ncid, id, self%psi_a)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
istatus = nf90_inq_varid (ncid, 'B_T', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='B_T')
istatus = nf90_get_var (ncid, id, self%B_T)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
! netcdf read scalars: R_mag,Z_mag,aminor
!
! R_mag == R position of magnetic axis
! Z_mag == Z position of magnetic axis
! aminor == half diameter of last flux surface at elevation of Z_mag
istatus = nf90_inq_varid (ncid, 'Rmag', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='Rmag')
istatus = nf90_get_var (ncid, id, self%R_mag)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
istatus = nf90_inq_varid (ncid, 'Zmag', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='Zmag')
istatus = nf90_get_var (ncid, id, self%Z_mag)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
istatus = nf90_inq_varid (ncid, 'a', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='a')
istatus = nf90_get_var (ncid, id, self%aminor)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
! netcdf read vectors: rho_d,eqpsi,psi_bar,fp,qsf,beta,pressure
!
! rho_d(1:nr) == half diameters of flux surfaces at elevation
! of Z_mag on the radial grid
! psi is the poloidal flux
! eqpsi(1:nr) == values of psi on the radial grid
! psi_bar(1:nr) == values of psi_bar on the radial grid
! [psi_bar == (eqpsi - psi_0)/(psi_a - psi_0) if not available]
! fp(1:nr) == the function that satisfies
! B = fp grad zeta + grad zeta x grad psi
! qsf(1:nr) == q profile on the radial grid
! beta(1:nr) == local beta, with the magnetic field defined
! to be vacuum magnetic field on axis
! pressure(1:nr) == pressure profile on the radial grid,
! normalized to the value at the magnetic axis.
istatus = nf90_inq_varid (ncid, 'psi', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='psi')
istatus = nf90_get_var (ncid, id, self%eqpsi)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
istatus = nf90_inq_varid (ncid, 'psibar', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='psibar')
istatus = nf90_get_var (ncid, id, self%psi_bar)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
istatus = nf90_inq_varid (ncid, 'fp', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='fp')
istatus = nf90_get_var (ncid, id, self%fp)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
istatus = nf90_inq_varid (ncid, 'q', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='q')
istatus = nf90_get_var (ncid, id, self%qsf)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
istatus = nf90_inq_varid (ncid, 'beta', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='beta')
istatus = nf90_get_var (ncid, id, self%beta)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
istatus = nf90_inq_varid (ncid, 'pressure', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='pressure')
istatus = nf90_get_var (ncid, id, self%pressure)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
! netcdf read 2d field: R_psi,Z_psi and B_psi (mod(B))
! eq_Rpsi(1:nr, 1:nt)
! eq_Zpsi(1:nr, 1:nt)
! eq_Bpsi(1:nr, 1:nt)
!
istatus = nf90_inq_varid (ncid, 'R_psi', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='R_psi')
istatus = nf90_get_var (ncid, id, self%R_psi, count = [self%nr, self%nt])
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
istatus = nf90_inq_varid (ncid, 'Z_psi', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='Z_psi')
istatus = nf90_get_var (ncid, id, self%Z_psi, count = [self%nr, self%nt])
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
istatus = nf90_inq_varid (ncid, 'B_psi', id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='B_psi')
istatus = nf90_get_var (ncid, id, self%B_psi, count = [self%nr, self%nt])
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)
istatus = nf90_close (ncid)
if (istatus /= NF90_NOERR) call netcdf_error (istatus)
! Convert units
self%psi_0 = self%psi_0 * psi_scale
self%psi_a = self%psi_a * psi_scale
self%eqpsi = self%eqpsi * psi_scale
self%B_T = self%B_T * b_scale ! cgs to MKS
self%B_psi = self%B_psi * b_scale / self%B_T
self%R_psi = self%R_psi / self%aminor
self%Z_psi = self%Z_psi / self%aminor
self%R_mag = self%R_mag / self%aminor
self%Z_mag = self%Z_mag / self%aminor
self%aminor = self%aminor * l_scale
self%fp = self%fp * b_scale * l_scale
self%beta_0 = self%beta(1)
! Normalize, rename quantities
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 = self%eqpsi / psi_N
f_N = self%B_T * self%aminor
self%fp = self%fp / f_N
do j = 1, self%nt
self%eqth(:, j) = (j-1) * pi / real(self%nt-1)
self%eqpsi_2d(:, j) = self%eqpsi
end do
self%diam = self%R_psi(:, 1) - self%R_psi(:, self%nt)
self%rc = 0.5 * (self%R_psi(:, 1) + self%R_psi(:, self%nt))
self%has_full_theta_range = .false.
# else
UNUSED_DUMMY(self) ; UNUSED_DUMMY(inputs)
call mp_abort('error: geq eqin is called without netcdf',.true.)
#endif
end subroutine read_and_set