submodule (polars_m) polars_s implicit none contains module subroutine set_parameters_polars(self, filename, nrho_in, ntheta_in, intorder_in) !! Sets parameters for polars, either via namelist from file, or setting explicitly implicit none class(polars_t), intent(inout) :: self character(*), intent(in), optional :: filename integer, intent(in), optional :: nrho_in integer, intent(in), optional :: ntheta_in integer, intent(in), optional :: intorder_in integer :: nrho, ntheta, intorder integer :: io_error character(len=256) :: io_errmsg namelist / polars / nrho, ntheta, intorder if ( (.not.present(filename)) .and. & (present(nrho_in)) .and. & (present(ntheta_in)) .and. & (present(intorder_in)) ) then ! input via explicit setting self%iwrk1 = nrho_in self%iwrk2 = ntheta_in self%intorder = intorder_in elseif ( (present(filename)) .and. & (.not.present(nrho_in)) .and. & (.not.present(ntheta_in)) .and. & (.not.present(intorder_in)) ) then ! Input via file ! default values intorder = 3 nrho = 0 ntheta = 0 open(unit = 20, file = filename, status = 'old', action = 'read',& iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif read(20, nml = polars, iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif close(20) self%iwrk1 = nrho self%iwrk2 = ntheta self%intorder = intorder else call handle_error('Input to subroutine not consistent', & GRILLIX_ERR_OTHER, __LINE__, __FILE__) endif end subroutine module subroutine build_polars(self, equi, mesh, dbgout) class(polars_t), intent(inout) :: self class(equilibrium_t) :: equi type(mesh_cart_t), intent(in) :: mesh integer, intent(in), optional :: dbgout integer :: nrho, ntheta, outi, ip, jp real(GP) :: phi, rho, theta phi = mesh%get_phi() ! set output level outi = 0 if (present(dbgout)) then if (is_rank_info_writer) then outi = dbgout endif if (dbgout >= 3) then outi = dbgout ! every rank writes endif endif nrho = self%iwrk1 ntheta = self%iwrk2 ! create polar grid call self%grid%initialize(equi, phi, nrho, ntheta) ! compute cartesian locations of polar grid allocate(self%xpol(ntheta, nrho)) allocate(self%ypol(ntheta, nrho)) !$omp parallel default(none) & !$omp private(ip, jp, rho, theta) & !$omp shared(self, equi, phi, nrho, ntheta) do ip = 1, nrho rho = self%grid%get_rho(ip) !$omp do do jp = 1, ntheta theta = self%grid%get_theta(jp) call polar_to_cart(equi, rho, theta, phi, self%xpol(jp, ip), self%ypol(jp, ip)) enddo !$omp end do enddo !$omp end parallel ! create map matrix cart -> polar call create_polar_map_matrix( self%cart_to_polar_csr, & equi, mesh, self%grid, self%intorder, outi) ! create surface average matrix call create_surface_average_csr( self%surface_average_csr, & equi, mesh, self%grid, self%cart_to_polar_csr, outi) ! create flux surface average matrix call create_flux_surface_average_csr( self%flux_surface_average_csr, & equi, mesh, self%grid, self%cart_to_polar_csr, outi) if (outi >0) then call self%display() endif ! Precompute fluxsurf areas and volumes if (.not. allocated(self%fluxsurf_area)) allocate(self%fluxsurf_area(nrho)) if (.not. allocated(self%fluxsurf_vol)) allocate(self%fluxsurf_vol(nrho)) do ip = 1, nrho self%fluxsurf_area(ip) = self%grid%fluxsurf_area(equi, ip) self%fluxsurf_vol(ip) = self%grid%fluxsurf_vol(equi, ip) end do ! nullify work variables self%iwrk1 = 0 self%iwrk2 = 0 end subroutine module subroutine display_polars(self) class(polars_t), intent(in) :: self if (.not. is_rank_info_writer) then return endif write(get_stdout(),*)'Information on polars ---------------------------------' call self%grid%display() write(get_stdout(),101) self%intorder 101 FORMAT( 3X,'intorder = ',I3 ) write(get_stdout(),*)'-------------------------------------------------------' end subroutine module subroutine write_netcdf_polars(self, equi, fgid) class(polars_t), intent(in) :: self class(equilibrium_t), intent(inout) :: equi integer, intent(in) :: fgid integer :: nf90_stat, id_grid, id_polmap, id_surfav, id_fluxsurfav, & iddim_nrho, id_fluxsurfarea, id_fluxsurfvol ! write global attributes (metadata) nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'intorder', self%intorder) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! Define dimensions nf90_stat = nf90_def_dim(fgid, 'nrho', self%grid%get_nrho(), & iddim_nrho) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! Define groups for grid and matrices nf90_stat = nf90_def_grp(fgid, 'grid', id_grid) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_def_grp(fgid, 'cart_to_polar', id_polmap) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_def_grp(fgid, 'surface_average', id_surfav) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_def_grp(fgid, 'flux_surface_average', id_fluxsurfav) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! Define variables nf90_stat = nf90_def_var(fgid, 'fluxsurf_area', NF90_DOUBLE, iddim_nrho, & id_fluxsurfarea) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_def_var(fgid, 'fluxsurf_vol', NF90_DOUBLE, iddim_nrho, & id_fluxsurfvol) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! End of definition nf90_stat = nf90_enddef(fgid) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! Write information call self%grid%write_netcdf(equi,id_grid) call self%cart_to_polar_csr%write_netcdf(id_polmap) call self%surface_average_csr%write_netcdf(id_surfav) call self%flux_surface_average_csr%write_netcdf(id_fluxsurfav) nf90_stat = nf90_put_var(fgid, id_fluxsurfarea, self%fluxsurf_area) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_put_var(fgid, id_fluxsurfvol, self%fluxsurf_vol) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) end subroutine module subroutine read_netcdf_polars(self, fgid) class(polars_t), intent(inout) :: self integer, intent(in) :: fgid integer :: nf90_stat, id_grid, id_polmap, id_surfav, id_fluxsurfav, & iddim_nrho, id_fluxsurfarea, id_fluxsurfvol, id_xpol, id_ypol real(GP), allocatable, dimension(:,:) :: wrk ! read global attributes nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'intorder', self%intorder) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! Read dimensions nf90_stat = nf90_inq_dimid(fgid, 'nrho', iddim_nrho) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! allocate matrices allocate(self%cart_to_polar_csr) allocate(self%surface_average_csr) allocate(self%flux_surface_average_csr) ! get groups for grid and matrices nf90_stat = nf90_inq_grp_ncid(fgid, 'grid', id_grid) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_inq_grp_ncid(fgid, 'cart_to_polar', id_polmap) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_inq_grp_ncid(fgid, 'surface_average', id_surfav) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_inq_grp_ncid(fgid, 'flux_surface_average', id_fluxsurfav) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! read grid and matrices call self%grid%read_netcdf(id_grid) call self%cart_to_polar_csr%read_netcdf(id_polmap) call self%surface_average_csr%read_netcdf(id_surfav) call self%flux_surface_average_csr%read_netcdf(id_fluxsurfav) ! Allocate fields allocate(self%fluxsurf_area(self%grid%get_nrho())) allocate(self%fluxsurf_vol(self%grid%get_nrho())) ! Read fields nf90_stat = nf90_inq_varid(fgid, 'fluxsurf_area', id_fluxsurfarea) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_get_var(fgid, id_fluxsurfarea, self%fluxsurf_area) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_inq_varid(fgid, 'fluxsurf_vol', id_fluxsurfvol) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_get_var(fgid, id_fluxsurfvol, self%fluxsurf_vol) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! Read Cartesian coordinates of polar grid ! these are inside the id_grid group, but in transposed ordering allocate(self%xpol(self%grid%get_ntheta(), self%grid%get_nrho())) allocate(self%ypol(self%grid%get_ntheta(), self%grid%get_nrho())) allocate(wrk(self%grid%get_nrho(), self%grid%get_ntheta())) nf90_stat = nf90_inq_varid(id_grid, 'xpol', id_xpol) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_get_var(id_grid, id_xpol, wrk) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) self%xpol = transpose(wrk) nf90_stat = nf90_inq_varid(id_grid, 'ypol', id_ypol) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_get_var(id_grid, id_ypol, wrk) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) self%ypol = transpose(wrk) end subroutine module subroutine destructor(self) type(polars_t), intent(inout) :: self if (allocated(self%xpol)) deallocate(self%xpol) if (allocated(self%ypol)) deallocate(self%ypol) if (allocated(self%cart_to_polar_csr)) deallocate(self%cart_to_polar_csr) if (allocated(self%surface_average_csr)) deallocate(self%surface_average_csr) if (allocated(self%flux_surface_average_csr)) deallocate(self%flux_surface_average_csr) if (allocated(self%fluxsurf_area)) deallocate(self%fluxsurf_area) if (allocated(self%fluxsurf_vol)) deallocate(self%fluxsurf_vol) self%iwrk1 = 0 self%iwrk2 = 0 self%intorder = 0 end subroutine end submodule