submodule (numerical_equilibrium_m) initialise_numerical_equilibrium_s use netcdf use bspline_module use bspline_kinds_module, only: wp use polygon_m use parallax_build_info_m, only: git_repo_path contains module subroutine read_params_in(self, filename) class(numerical_t), intent(inout) :: self character(len = *), intent(in), OPTIONAL :: filename integer :: io_error character(len = string_buffer_length) :: equilibrium_case character(len = string_buffer_length) :: path_to_netcdf logical :: flip_Z integer :: spline_order real(FP) :: rhomin, rhomax, rhomin_privflux, poloidal_field_factor logical :: set_rhomin, set_rhomax, set_rhomin_privflux, raise_exception_off_grid character(len=:), allocatable :: district_definition logical :: set_exclusion_polygon, set_divertor_polygon, invert_exclusion_polygon, invert_divertor_polygon ! Buffer for manually setting polygons integer, parameter :: max_polygon_N_pts = 100 integer :: divertor_polygon_N_pts, exclusion_polygon_N_pts real(FP), dimension(max_polygon_N_pts) :: divertor_polygon_X_pts real(FP), dimension(max_polygon_N_pts) :: divertor_polygon_Y_pts real(FP), dimension(max_polygon_N_pts) :: exclusion_polygon_X_pts real(FP), dimension(max_polygon_N_pts) :: exclusion_polygon_Y_pts integer :: iunit = 35 character(len = string_buffer_length), allocatable :: line namelist / equi_numerical_params /equilibrium_case , & & path_to_netcdf , & & flip_Z , & & poloidal_field_factor , & & spline_order , & & rhomin, rhomax, rhomin_privflux, & & raise_exception_off_grid, & & district_definition, & & set_rhomin, set_rhomax, set_rhomin_privflux, & & set_exclusion_polygon, set_divertor_polygon, & & invert_exclusion_polygon, invert_divertor_polygon, & divertor_polygon_N_pts, divertor_polygon_X_pts, divertor_polygon_Y_pts, & exclusion_polygon_N_pts, exclusion_polygon_X_pts, exclusion_polygon_Y_pts ! Make some friendly default values spline_order = 4 set_rhomin = .false. set_rhomax = .false. set_rhomin_privflux = .false. raise_exception_off_grid = .false. district_definition = 'flux' rhomin = 0.0_FP rhomax = 0.0_FP rhomin_privflux = 0.0_FP set_divertor_polygon = .false. set_exclusion_polygon = .false. invert_exclusion_polygon = .false. invert_divertor_polygon = .false. flip_Z = .false. poloidal_field_factor = 1.0_FP divertor_polygon_N_pts = 0 divertor_polygon_X_pts = 0.0_FP divertor_polygon_Y_pts = 0.0_FP exclusion_polygon_N_pts = 0 exclusion_polygon_X_pts = 0.0_FP exclusion_polygon_Y_pts = 0.0_FP if (present(filename)) then open(unit = iunit, file = filename, status = 'old', action = 'read', iostat = io_error) if (io_error /= 0) then write(get_stderr(), *) & 'error(equi%init): file', filename, 'not existent' ERROR STOP ERR_PARAMETER_FILE endif read(iunit, nml = equi_numerical_params, iostat = io_error) if (io_error == -1) then write(get_stderr(),'(A)') & 'End-of-file reached: either missing a parameter group, & &or missing a blank line at the end of file' write(get_stderr(), *) & 'error(equi%init): erroneous namelist', io_error ERROR STOP ERR_PARAMETER_FILE elseif(io_error /= 0) then backspace(iunit) read(iunit,fmt='(A)') line write(get_stderr(),'(A)') & 'Invalid line in namelist: '//trim(line) write(get_stderr(), *) & 'error(equi%init): erroneous namelist', io_error ERROR STOP ERR_PARAMETER_FILE endif close(iunit) else ! If not given a parameter file, set default for CI/CD equilibrium_case = "SF_fourier_sigma=2d00_original" path_to_netcdf = git_repo_path // "/res/numerical_equilibria" set_rhomin = .true. rhomin = 0.90 set_rhomax = .true. rhomax = 1.02 endif ! Copy params into the equi object self%equilibrium_case = equilibrium_case self%path_to_netcdf = path_to_netcdf self%flip_Z = flip_Z self%poloidal_field_factor = poloidal_field_factor self%spline_order = spline_order self%district_definition = district_definition self%set_rhomin = set_rhomin self%set_rhomax = set_rhomax self%set_rhomin_privflux = set_rhomin_privflux self%rhomin = rhomin self%rhomax = rhomax self%rhomin_privflux_params = rhomin_privflux self%raise_exception_off_grid = raise_exception_off_grid self%set_divertor_polygon = set_divertor_polygon self%set_exclusion_polygon = set_exclusion_polygon self%invert_divertor_polygon = invert_divertor_polygon self%invert_exclusion_polygon = invert_exclusion_polygon if (self%set_divertor_polygon) then if (self%flip_Z) then ! Reverse order to keep anti-clockwise, Z -> -Z call self%make_polygon(N_points = divertor_polygon_N_pts, & X_points = divertor_polygon_X_pts(divertor_polygon_N_pts:1:-1), & Y_points = -1 * divertor_polygon_Y_pts(divertor_polygon_N_pts:1:-1), & polygon = self%divertor_polygon) else call self%make_polygon(N_points = divertor_polygon_N_pts, & X_points = divertor_polygon_X_pts(1:divertor_polygon_N_pts), & Y_points = divertor_polygon_Y_pts(1:divertor_polygon_N_pts), & polygon = self%divertor_polygon) endif endif if (self%set_exclusion_polygon) then if (self%flip_Z) then ! Reverse order to keep anti-clockwise, Z -> -Z call self%make_polygon(N_points = exclusion_polygon_N_pts, & X_points = exclusion_polygon_X_pts(exclusion_polygon_N_pts:1:-1), & Y_points = -1 * exclusion_polygon_Y_pts(exclusion_polygon_N_pts:1:-1), & polygon = self%exclusion_polygon) else call self%make_polygon(N_points = exclusion_polygon_N_pts, & X_points = exclusion_polygon_X_pts(1:exclusion_polygon_N_pts), & Y_points = exclusion_polygon_Y_pts(1:exclusion_polygon_N_pts), & polygon = self%exclusion_polygon) endif endif self%params_filename = filename end subroutine read_params_in module subroutine read_from_ncdf(self) class(numerical_t), intent(inout) :: self integer :: root_grp, status character(len = string_buffer_length) :: full_filepath integer :: geo_grp, psi_lim_grp integer :: grid_lims_grp, divertor_polygon_grp, exclusion_polygon_grp write(full_filepath, '(A)') trim(self%path_to_netcdf) // "/" // trim(self%equilibrium_case) // ".nc" !$OMP CRITICAL ! Load a NetCDF located in the run-time directory status = nf90_open(path = trim(full_filepath), mode = nf90_nowrite, ncid = root_grp) if (status /= nf90_noerr) then write(get_stderr(), *) & "Error. Failed to open equilibrium netCDF. & &Filepath given was ", trim(full_filepath) ERROR STOP ERR_UNHANDLED endif if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_from_ncdf::nf90_open") status = nf90_inq_ncid(root_grp, "Magnetic_geometry", geo_grp) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_from_ncdf::nf90_inq_ncid for geo_grp") status = nf90_inq_ncid(root_grp, "Psi_limits", psi_lim_grp) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_from_ncdf::nf90_inq_ncid for psi_lim_grp") status = nf90_inq_ncid(root_grp, "Grid_limits", grid_lims_grp) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_from_ncdf::nf90_inq_ncid for grid_lims_grp") status = nf90_inq_ncid(root_grp, "divertor_polygon", divertor_polygon_grp) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_from_ncdf::nf90_inq_ncid for divertor_polygon_grp") status = nf90_inq_ncid(root_grp, "exclusion_polygon", exclusion_polygon_grp) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_from_ncdf::nf90_inq_ncid for exclusion_polygon_grp") call find_ncdf_attributes(self, root_grp) call make_interpolators(self, geo_grp) call find_psi_limits(self, psi_lim_grp) call find_box_limits(self, grid_lims_grp) if (.not.(self%set_divertor_polygon)) then ! If set_divertor_polygon parameter is .true., it should be set in the ! read parameters step. Otherwise (default option) read the polygon ! from the NetCDF call extract_nc_polygon(self, self%divertor_polygon, divertor_polygon_grp) endif if (.not.(self%set_exclusion_polygon)) then ! If set_exclusion_polygon parameter is .true., it should be set in the ! read parameters step. Otherwise (default option) read the polygon ! from the NetCDF call extract_nc_polygon(self, self%exclusion_polygon, exclusion_polygon_grp) endif status = nf90_close(ncid = root_grp) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_from_ncdf::nf90_close") !$OMP END CRITICAL end subroutine subroutine find_ncdf_attributes(self, root_grp) class(numerical_t), intent(inout) :: self integer, intent(in) :: root_grp integer :: status status = nf90_get_att(root_grp, nf90_global, "description", self%ncdf_description) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_ncdf_attributes::nf90_get_att") status = nf90_get_att(root_grp, nf90_global, "history", self%ncdf_history) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_ncdf_attributes::nf90_get_att") status = nf90_get_att(root_grp, nf90_global, "author", self%ncdf_author) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_ncdf_attributes::nf90_get_att") end subroutine find_ncdf_attributes subroutine make_interpolators(self, geo_grp) class(numerical_t), intent(inout) :: self integer, intent(in) :: geo_grp integer :: Nr, Nz real(FP), dimension(:), allocatable :: R, Z real(FP), dimension(5) :: X_points, Y_points call interrogate_axis_in_group(geo_grp, "R", Nr, R) call interrogate_axis_in_group(geo_grp, "Z", Nz, Z) call read_critical_points(self, geo_grp) ! Normalise the axes to the magnetic axis R R = R if (self%flip_Z) then ! Reverse the order of the Z values so that they are still ascending order Z = -1 * Z(Nz:1:-1) else Z = Z endif call write_spline_info(self, R, Z, Nr, Nz) ! Build spline limit points (without literals) X_points(1) = self%spl_limits_R(2) Y_points(1) = self%spl_limits_Z(2) X_points(2) = self%spl_limits_R(1) Y_points(2) = self%spl_limits_Z(2) X_points(3) = self%spl_limits_R(1) Y_points(3) = self%spl_limits_Z(1) X_points(4) = self%spl_limits_R(2) Y_points(4) = self%spl_limits_Z(1) X_points(5) = self%spl_limits_R(2) Y_points(5) = self%spl_limits_Z(2) ! Have already flipped Z array, order is defined in terms of min, max -- reversal is automatic call self%make_polygon(5, X_points, Y_points, self%spline_edge_polygon) ! Have already flipped Z array call build_spline_interpolant(geo_grp, "psi", self%psi_spl2d, R, Z, Nr, Nz, self%spline_order, self%flip_Z) deallocate(R) deallocate(Z) end subroutine make_interpolators subroutine interrogate_axis_in_group(ncid, name, dim_length, dim_values) integer, intent(in) :: ncid character(len=*), intent(in) :: name integer, intent(out) :: dim_length real(FP), allocatable, dimension(:), intent(out) :: dim_values integer :: dimid, varid, status ! Find the dimension ID status = nf90_inq_dimid(ncid, name, dimid) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &interrogate_axis::nf90_inq_dimid") ! Find the length of the dimension status = nf90_inquire_dimension(ncid, dimid, len = dim_length) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &interrogate_axis::nf90_inquire_dimension") ! Find the variable ID for the dimension values status = nf90_inq_varid(ncid, name, varid) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &interrogate_axis::nf90_inq_varid") ! Allocate memory for the dimension values allocate(dim_values(dim_length)) ! Fill the dimension values status = nf90_get_var(ncid, varid, values = dim_values) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &interrogate_axis::nf90_get_var") end subroutine subroutine read_critical_points(self, geo_grp) ! Normalise all spatial values to R0 -- this is done in preprocessing, except for R0 and Z0 class(numerical_t), intent(inout) :: self integer, intent(in) :: geo_grp integer :: status real(FP) :: R0, Z0 status = nf90_get_att(geo_grp, nf90_global, "magnetic_axis_R", R0) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_critical_points::nf90_get_att magnetic_axis_R") self%R0 = R0 self%x0 = 1.0_FP status = nf90_get_att(geo_grp, nf90_global, "magnetic_axis_Z", Z0) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_critical_points::nf90_get_att magnetic_axis_Z") if (self%flip_Z) then self%y0 = -1 * Z0/self%R0 else self%y0 = Z0/self%R0 endif status = nf90_get_att(geo_grp, nf90_global, "axis_Btor", self%axis_Btor) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_critical_points::nf90_get_att axis_Btor") ! if (self%axis_Btor < 0 .and. .not.(self%flip_Z)) then ! if(is_master()) write(stdout, *) "Warning (initialise_numerical_equilibrium_m)::& ! &for consistency with experiment, recommend & ! &using flip_Z = .true. to reverse the toroidal magnetic field. & ! &self%axis_Btor is taken always >0." ! endif self%axis_Btor = abs(self%axis_Btor) end subroutine read_critical_points subroutine build_spline_interpolant(ncid, name, spline2D, R, Z, Nr, Nz, spline_order, flip_Z) integer, intent(in) :: ncid character(len=*), intent(in) :: name type(bspline_2d), intent(inout) :: spline2D integer :: varid, status integer, intent(in) :: Nr, Nz real(FP), dimension(Nr), intent(in) :: R real(FP), dimension(Nz), intent(in) :: Z real(FP), dimension(Nr, Nz) :: nc_values integer, intent(in) :: spline_order logical, intent(in) :: flip_Z CHARACTER(len=:), ALLOCATABLE :: msg ! Get variable ID from the NetCDF status = nf90_inq_varid(ncid, name, varid) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &build_spline_interpolant::nf90_inq_varid") status = nf90_get_var(ncid, varid, values = nc_values) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &build_spline_interpolant::nf90_get_var") if (flip_Z) then ! Need to flip the field values call spline2D%initialize(R, Z, nc_values(1:Nr:1, Nz:1:-1), spline_order, spline_order, status) else call spline2D%initialize(R, Z, nc_values, spline_order, spline_order, status) endif ! Uncomment to allow extrapolation (not recommended) ! call spline2D%set_extrap_flag(.true.) if (status /= 0) then call get_status_message(status, msg) write(get_stderr(),*) 'Error initializing 2D spline: '//msg ERROR STOP ERR_UNHANDLED endif end subroutine subroutine write_spline_info(self, R, Z, Nr, Nz) class(numerical_t), intent(inout) :: self integer, intent(in) :: Nr, Nz real(FP), dimension(Nr), intent(in) :: R real(FP), dimension(Nz), intent(in) :: Z ! Write information about the spline interpolator to module variables ! N.b. R and Z are normalised to the major radius self%spl_limits_R(1) = minval(R) self%spl_limits_R(2) = maxval(R) self%spl_size_R = Nr self%spl_limits_Z(1) = minval(Z) self%spl_limits_Z(2) = maxval(Z) self%spl_size_Z = Nz end subroutine write_spline_info subroutine find_psi_limits(self, psi_lim_grp) class(numerical_t), intent(inout) :: self integer, intent(in) :: psi_lim_grp call find_global_limits(self, psi_lim_grp) call find_region_limits(self, psi_lim_grp) end subroutine find_psi_limits subroutine find_global_limits(self, psi_lim_grp) class(numerical_t), intent(inout) :: self integer, intent(in) :: psi_lim_grp integer :: status integer :: xtype, len, attnum ! We will only use len (length). xtype gives the type of the output, and attnum is the attribute number real(FP), ALLOCATABLE, DIMENSION(:) :: seperatrices_psi if (.not.(self%set_rhomin)) then status = nf90_get_att(psi_lim_grp, nf90_global, "rho_min", self%rhomin) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_global_limits::nf90_get_att rhomin") endif if (.not.(self%set_rhomax)) then status = nf90_get_att(psi_lim_grp, nf90_global, "rho_max", self%rhomax) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_global_limits::nf90_get_att rhomax") endif ! It is possible to have multiple seperatrices defined in the NetCDF. We want the maximum (innermost) status = nf90_inquire_attribute(psi_lim_grp, nf90_global, "psi_seperatrix", xtype, len, attnum) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_critical_points::nf90_inquire_attribute psi_seperatrix") ALLOCATE(seperatrices_psi(1:len)) status = nf90_get_att(psi_lim_grp, nf90_global, "psi_seperatrix", seperatrices_psi) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_critical_points::nf90_get_att psi_seperatrix") self%X_point_psi = maxval(seperatrices_psi) DEALLOCATE(seperatrices_psi) status = nf90_get_att(psi_lim_grp, nf90_global, "psi_axis", self%O_point_psi) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &read_critical_points::nf90_get_att psi_axis") end subroutine find_global_limits subroutine find_region_limits(self, psi_lim_grp) class(numerical_t), intent(inout) :: self integer, intent(in) :: psi_lim_grp integer :: id_regions(flux_regions_buffer) integer :: status, region_it real(FP) :: local_min, local_max integer :: use_local_min_int, use_local_max_int logical :: use_local_min, use_local_max id_regions = 0 status = nf90_inq_grps(psi_lim_grp, self%num_flux_regions, id_regions) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_region_limits::nf90_inq_grps") allocate(self%flux_regions(self%num_flux_regions)) do region_it = 1, self%num_flux_regions status = nf90_get_att(id_regions(region_it), nf90_global, "local_rho_min", local_min) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_region_limits::nf90_get_att local_min") status = nf90_get_att(id_regions(region_it), nf90_global, "local_rho_max", local_max) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_region_limits::nf90_get_att local_max") status = nf90_get_att(id_regions(region_it), nf90_global, "use_local_min", use_local_min_int) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_region_limits::nf90_get_att use_local_min") status = nf90_get_att(id_regions(region_it), nf90_global, "use_local_max", use_local_max_int) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_region_limits::nf90_get_att use_local_max") ! Make sure that the integers read correctly if (use_local_min_int == 1) then use_local_min = .true. if (self%set_rhomin_privflux) then local_min = self%rhomin_privflux_params endif elseif (use_local_min_int == 0) then use_local_min = .false. else write(get_stderr(), *) & "Error(initialise_numerical_equilibrium_t::& &find_region_limits). use_local_min = ", use_local_min_int endif if (use_local_max_int == 1) then use_local_max = .true. elseif (use_local_max_int == 0) then use_local_max = .false. else write(get_stderr(), *) & "Error(initialise_numerical_equilibrium_t::& &find_region_limits). use_local_max = ", use_local_max_int endif call extract_nc_polygon(self, self%flux_regions(region_it), id_regions(region_it), & local_min, local_max, use_local_min, use_local_max) end do end subroutine find_region_limits subroutine find_box_limits(self, grid_lims_grp) class(numerical_t), intent(inout) :: self integer, intent(in) :: grid_lims_grp integer :: status real(FP) :: rmin, rmax, zmin, zmax status = nf90_get_att(grid_lims_grp, nf90_global, "grid_rmin", rmin) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_global_limits::nf90_get_att") status = nf90_get_att(grid_lims_grp, nf90_global, "grid_rmax", rmax) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_global_limits::nf90_get_att") status = nf90_get_att(grid_lims_grp, nf90_global, "grid_zmin", zmin) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_global_limits::nf90_get_att") status = nf90_get_att(grid_lims_grp, nf90_global, "grid_zmax", zmax) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &find_global_limits::nf90_get_att") ! Check if undefined self%xmin = rmin self%xmax = rmax self%ymin = zmin self%ymax = zmax if (self%flip_Z) then zmin = self%ymin zmax = self%ymax self%ymin = -zmax self%ymax = -zmin endif end subroutine find_box_limits subroutine extract_nc_polygon(self, nc_polygon, grp, local_min, local_max, use_local_min, use_local_max) class(numerical_t), intent(inout) :: self integer, intent(in) :: grp real(FP), intent(in), optional :: local_min, local_max logical, intent(in), optional :: use_local_min, use_local_max class(polygon2d_t), intent(out) :: nc_polygon real(FP), dimension(:), allocatable :: R_points, Z_points integer :: N_points call find_RZ_points_from_nc(grp, R_points, Z_points, N_points) if (self%flip_Z) then R_points = R_points(N_points:1:-1) Z_points = -1 * Z_points(N_points:1:-1) endif select type (nc_polygon) type is (closed_polygon2d_t) nc_polygon = create_closed_polygon2d_t(N_points, R_points, Z_points) type is (limiting_polygon2d_t) if(present(local_min) .and. present(local_max) .and. present(use_local_max) .and. present(use_local_max)) then nc_polygon = create_limiting_polygon2d_t(N_points, R_points, Z_points, & local_min, local_max, use_local_min, use_local_max) else write(get_stderr(), *) & "error: limiting_polygon2d_t called without & &defined limits" ERROR STOP ERR_UNHANDLED endif class default write(get_stderr(), *) "error: nc_polygon class not recognised" ERROR STOP ERR_UNHANDLED end select deallocate(R_points, Z_points) end subroutine extract_nc_polygon subroutine find_RZ_points_from_nc(grp, R_points, Z_points, dim_length) integer, intent(in) :: grp ! Has dimension N_points ! Variables R_points, Z_points real(FP), dimension(:), allocatable, intent(inout) :: R_points, Z_points integer :: dimid, dim_length integer :: Rid, Zid, status ! Find the dimension ID status = nf90_inq_dimid(grp, "N_points", dimid) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &extract_nc_polygon::nf90_inq_dimid") ! Find the length of the dimension status = nf90_inquire_dimension(grp, dimid, len = dim_length) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &extract_nc_polygon::nf90_inquire_dimension") ! Find the variable ID for the polygon points status = nf90_inq_varid(grp, "R_points", Rid) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &extract_nc_polygon::nf90_inq_varid") status = nf90_inq_varid(grp, "Z_points", Zid) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &extract_nc_polygon::nf90_inq_varid") ! Allocate memory for the polygon points allocate(R_points(dim_length)) allocate(Z_points(dim_length)) ! Fill the polygon points status = nf90_get_var(grp, Rid, values = R_points) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &extract_nc_polygon::nf90_get_var") status = nf90_get_var(grp, Zid, values = Z_points) if (status /= nf90_noerr) call nf90_handle_err(status, "numerical_equilibrium_t::& &extract_nc_polygon::nf90_get_var") end subroutine find_RZ_points_from_nc end submodule initialise_numerical_equilibrium_s