read_and_set_teq Subroutine

  subroutine read_and_set_teq(self, inputs)
    implicit none
    class(teq_type), intent(in out) :: self
    type(geo_input_type), intent(in) :: inputs
# ifdef NETCDF
    integer :: istatus, ncid, id
    real, allocatable, dimension(:,:) :: work
    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)

    !     nr == number of radial grid points in radial eq grid
    istatus = nf90_inq_varid (ncid, 'ns', id)
    if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='ns')
    istatus = nf90_get_var (ncid, id, self%nr)
    if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)

    !     nt == number of theta grid points in theta eq grid
    istatus = nf90_inq_varid (ncid, 'nt1', id)
    if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='nt1')
    istatus = nf90_get_var (ncid, id, self%nt)
    if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)

    call self%alloc_arrays()

    !     netcdf read vectors: eqpsi,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, 'g', id)
    if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='g')
    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, 'p', id)
    if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='p')
    istatus = nf90_get_var (ncid, id, self%pressure)
    if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)

    allocate(work(self%nt, self%nr))

    !     netcdf read 2d field: R_psi and Z_psi
    istatus = nf90_inq_varid (ncid, 'R', id)
    if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='R')
    istatus = nf90_get_var (ncid, id, work)
    if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)

    self%R_psi = transpose(work)

    istatus = nf90_inq_varid (ncid, 'Z', id)
    if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='Z')
    istatus = nf90_get_var (ncid, id, work)
    if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, id)

    self%Z_psi = transpose(work)

    istatus = nf90_close (ncid)
    if (istatus /= NF90_NOERR) call netcdf_error (istatus)
    deallocate(work)

    self%eqpsi = abs(self%eqpsi) ; self%fp = abs(self%fp) ; self%qsf = abs(self%qsf)
    self%has_full_theta_range = .true.
    call self%shared_setup(self%nt, self%nt/2)
# else
    UNUSED_DUMMY(self) ; UNUSED_DUMMY(inputs)
    call mp_abort('error: peq teqin is called without netcdf',.true.)
# endif
  end subroutine read_and_set_teq