submodule(dommaschk_equilibrium_m) dommaschk_equilibrium_netcdf_s use netcdf use array_generation_m, only : linspace use error_handling_m, only : handle_error_netcdf implicit none contains module subroutine read_bndry_polygons(self, filename, variable, & phi_array, polygon_vertices) class(dommaschk_t), intent(in) :: self character(len=*), intent(in) :: filename character(len=*), intent(in) :: variable real(FP), dimension(:), allocatable :: phi_array real(FP), dimension(:,:,:), allocatable, intent(out) :: polygon_vertices integer :: file_id, ierr integer :: id_dim_n_vertices, id_dim_n_planes, n_vertices, n_planes integer :: id_var_phi_array, id_var_polygon ierr = nf90_open(trim(filename), nf90_nowrite, file_id) call handle_error_netcdf(ierr, __LINE__, __FILE__) ! Dimensions ierr = nf90_inq_dimid(file_id, 'n_vertices', id_dim_n_vertices) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_inquire_dimension(file_id, id_dim_n_vertices, len=n_vertices) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_inq_dimid(file_id, 'n_planes', id_dim_n_planes) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_inquire_dimension(file_id, id_dim_n_planes, len=n_planes) call handle_error_netcdf(ierr, __LINE__, __FILE__) ! Read data allocate(phi_array(n_planes)) allocate(polygon_vertices(n_vertices, 2, n_planes)) ierr = nf90_inq_varid(file_id, 'phi_array', id_var_phi_array) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_get_var(file_id, id_var_phi_array, phi_array) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_inq_varid(file_id, trim(variable), id_var_polygon) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_get_var(file_id, id_var_polygon, polygon_vertices) call handle_error_netcdf(ierr, __LINE__, __FILE__) ! Close file ierr = nf90_close(file_id) call handle_error_netcdf(ierr, __LINE__, __FILE__) if (n_planes > 1) then call check_phi_array_consistency(self%num_field_periods, n_planes, & phi_array) endif ! Ensure polygons are closed if(maxval(abs( polygon_vertices(1, :, :) & - polygon_vertices(n_vertices, :, :))) > FP_EPS) then call handle_error("Polygon "//trim(variable)//" read from file "//& trim(filename)//" is not closed!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__) endif end subroutine read_bndry_polygons module subroutine read_rho_polygons(self, filename) class(dommaschk_t), intent(inout) :: self character(len=*), intent(in) :: filename integer :: file_id, ierr integer :: id_dim_n_vertices, id_dim_n_planes, id_dim_n_surf, & id_dim_n_surf_excluded integer :: n_surfaces, n_surf_excluded, n_vertices, n_planes integer :: id_var_phi_array, id_var_polygon, id_var_surf_excluded integer, dimension(:), allocatable :: surf_excluded integer :: i_surf, i_included character(len=256) :: variable ierr = nf90_open(trim(filename), nf90_nowrite, file_id) call handle_error_netcdf(ierr, __LINE__, __FILE__) ! Dimensions ierr = nf90_inq_dimid(file_id, 'n_vertices', id_dim_n_vertices) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_inquire_dimension(file_id, id_dim_n_vertices, len=n_vertices) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_inq_dimid(file_id, 'n_planes', id_dim_n_planes) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_inquire_dimension(file_id, id_dim_n_planes, len=n_planes) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_inq_dimid(file_id, 'n_surf', id_dim_n_surf) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_inquire_dimension(file_id, id_dim_n_surf, len=n_surfaces) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_inq_dimid(file_id, 'n_surf_excluded', id_dim_n_surf_excluded) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_inquire_dimension(file_id, id_dim_n_surf_excluded, & len=n_surf_excluded) call handle_error_netcdf(ierr, __LINE__, __FILE__) ! Read data allocate(self%rho_polygons_phi(n_planes)) allocate(surf_excluded(n_surf_excluded)) allocate(self%rho_polygons(n_vertices, 2, n_planes, n_surfaces)) ierr = nf90_inq_varid(file_id, 'phi_array', id_var_phi_array) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_get_var(file_id, id_var_phi_array, self%rho_polygons_phi) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_inq_varid(file_id, 'surf_excluded', id_var_surf_excluded) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_get_var(file_id, id_var_surf_excluded, surf_excluded) call handle_error_netcdf(ierr, __LINE__, __FILE__) i_included = 1 do i_surf = 0, n_surfaces + n_surf_excluded - 1 ! Some surfaces are excluded from the Poincare data because they do not ! form smooth polygons (rational surfaces, islands, etc.) if(any(i_surf == surf_excluded)) then cycle endif write(variable, "(A, I3.3)") "polygon_", i_surf ierr = nf90_inq_varid(file_id, trim(variable), id_var_polygon) call handle_error_netcdf(ierr, __LINE__, __FILE__) ierr = nf90_get_var(file_id, id_var_polygon, self%rho_polygons(:,:,:,i_included)) call handle_error_netcdf(ierr, __LINE__, __FILE__) i_included = i_included + 1 enddo ! Close file ierr = nf90_close(file_id) call handle_error_netcdf(ierr, __LINE__, __FILE__) if (n_planes > 1) then call check_phi_array_consistency(self%num_field_periods, n_planes, & self%rho_polygons_phi) endif ! Ensure polygons are closed if(maxval(abs( self%rho_polygons(1, :, :, :) & - self%rho_polygons(n_vertices, :, :, :))) > FP_EPS) then call handle_error("Flux surface polygons read from file "//& trim(filename)//" are not closed!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__) endif end subroutine read_rho_polygons subroutine check_phi_array_consistency(num_field_periods, n_planes, phi_array) !! Check that phi_array is evenly spaced in first field period and !! sorted in increasing order integer, intent(in) :: num_field_periods !! Number of field periods in equilibrium integer, intent(in) :: n_planes !! Number of toroidal planes real(FP), dimension(n_planes) :: phi_array !! Toroidal angles integer :: k real(FP), dimension(:), allocatable :: phi_array_expected allocate(phi_array_expected(n_planes)) phi_array_expected = linspace(0.0_FP, & TWO_PI / num_field_periods, & n_planes, endpoint=.false.) do k = 1, n_planes if (abs(phi_array(k) - phi_array_expected(k)) > FP_EPS) then call handle_error("Polygons are not defined at expected & &phi values!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__, & additional_info=& error_info_t("is, expected:", & [phi_array(k), phi_array_expected(k)])) endif enddo deallocate(phi_array_expected) end subroutine end submodule