submodule(polar_grid_m) polar_grid_s use screen_io_m, only : get_stderr, get_stdout, nf90_handle_err implicit none contains module subroutine initialize_polar_grid(self, equi, phi, nrho, ntheta) class(polar_grid_t), intent(inout) :: self class(equilibrium_t), intent(inout) :: equi real(FP), intent(in) :: phi integer, intent(in) :: nrho integer, intent(in) :: ntheta self%nrho = nrho self%ntheta = ntheta self%phi = phi if ((nrho == 0) .or. (ntheta == 0)) then ! For empty polar grid self%drho = FP_NAN self%dtheta = FP_NAN return endif select type(equi) type is(circular_t) self%rhopol_min = equi%rhomin self%rhopol_max = equi%rhomax type is (slab_t) self%rhopol_min = equi%rhomin self%rhopol_max = equi%rhomax class default ! Establish polar grid only for closed flux surface region if (equi%rhomin >= 1.0_FP) then self%rhopol_min = 0.0_FP self%rhopol_max = 0.0_FP self%ntheta = 0 self%nrho = 0 self%drho = 0.0_FP self%dtheta = 0.0_FP return else self%rhopol_min = equi%rhomin if (equi%rhomax < 1.0_FP) then self%rhopol_max = equi%rhomax else ! rhopolmax = 1-drho self%drho = (1.0_FP-equi%rhomin) / (1.0_FP*nrho) self%rhopol_max = self%rhopol_min + (nrho-1)*self%drho endif endif end select self%drho = (self%rhopol_max - self%rhopol_min) / (1.0_FP*(nrho-1)) self%dtheta = two_pi / (1.0*ntheta) end subroutine pure real(FP) module function get_phi(self) class(polar_grid_t), intent(in) :: self get_phi = self%phi end function pure real(FP) module function get_rhopol_min(self) class(polar_grid_t), intent(in) :: self get_rhopol_min = self%rhopol_min end function pure real(FP) module function get_rhopol_max(self) class(polar_grid_t), intent(in) :: self get_rhopol_max = self%rhopol_max end function pure integer module function get_nrho(self) class(polar_grid_t), intent(in) :: self get_nrho = self%nrho end function pure integer module function get_ntheta(self) class(polar_grid_t), intent(in) :: self get_ntheta = self%ntheta end function pure real(FP) module function get_drho(self) class(polar_grid_t), intent(in) :: self get_drho = self%drho end function pure real(FP) module function get_dtheta(self) class(polar_grid_t), intent(in) :: self get_dtheta = self%dtheta end function pure real(FP) module function get_rho(self, ip) class(polar_grid_t), intent(in) :: self integer, intent(in) :: ip real(FP), parameter :: singularity_corr_fac = 1.0E-4_FP get_rho = self%rhopol_min+(ip-1)*self%drho ! Circumvent singularity at magnetic axis if ( (ip == 1) .and. (get_rho < FP_EPS) ) then get_rho = get_rho + singularity_corr_fac * self%drho endif end function pure real(FP) module function get_theta(self, jp) class(polar_grid_t), intent(in) :: self integer, intent(in) :: jp get_theta = (jp-1)*self%dtheta end function real(FP) module function fluxsurf_area(self, equi, ip) class(polar_grid_t), intent(in) :: self class(equilibrium_t), intent(inout) :: equi integer, intent(in) :: ip integer :: jp real(FP) :: rho, theta, phi, dtheta, area real(FP) :: xp, yp, xl1, xl2, yl1, yl2, lpol, ltor phi = self%get_phi() dtheta = self%get_dtheta() area = 0.0_FP do jp = 1, self%ntheta rho = self%get_rho(ip) theta = self%get_theta(jp) call polar_to_cart(equi, rho, theta, phi, xp, yp) call polar_to_cart(equi, rho, theta - dtheta / 2.0_FP, phi, xl1, yl1) call polar_to_cart(equi, rho, theta + dtheta / 2.0_FP, phi, xl2, yl2) lpol = sqrt((xl1 - xl2)**2 + (yl1 - yl2)**2) ltor = TWO_PI * equi%jacobian(xp, yp, phi) area = area + lpol*ltor enddo fluxsurf_area = area end function real(FP) module function fluxsurf_vol(self, equi, ip) class(polar_grid_t), intent(in) :: self class(equilibrium_t), intent(inout) :: equi integer, intent(in) :: ip integer :: jp real(FP) :: rho, theta, phi, dtheta, drho, vol real(FP) :: xp, yp phi = self%get_phi() drho = self%get_drho() dtheta = self%get_dtheta() vol = 0.0_FP do jp = 1, self%ntheta rho = self%get_rho(ip) theta = self%get_theta(jp) call polar_to_cart(equi, rho, theta, phi, xp, yp) vol = vol + TWO_PI * drho * dtheta * & jacobian_polar(equi, xp, yp, phi) enddo fluxsurf_vol = vol end function module subroutine write_netcdf_polar(self, equi, fgid) class(polar_grid_t), intent(in) :: self class(equilibrium_t), intent(inout) :: equi integer, intent(in) :: fgid integer :: nf90_stat, iddim_nrho, iddim_ntheta, id_rho, id_theta, id_xpol, id_ypol integer :: ip, jp real(FP), dimension(self%nrho) :: rho_arr real(FP), dimension(self%ntheta) :: theta_arr real(FP), dimension(self%nrho, self%ntheta) :: xpol, ypol character(len=*), parameter :: cname = "write_netcdf_polar" ! Name used for logging errors ! Write global attributes (metadata) ------------------------------------------------------ nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'phi', self%phi) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'rhopol_min', self%rhopol_min) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'rhopol_max', self%rhopol_max) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'drho', self%drho) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'dtheta', self%dtheta) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) ! Define dimensions ----------------------------------------------------------------------- nf90_stat = nf90_def_dim(fgid, 'nrho', self%nrho, iddim_nrho) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_def_dim(fgid, 'ntheta', self%ntheta, iddim_ntheta) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) ! Define variables ------------------------------------------------------------------------ nf90_stat = nf90_def_var(fgid, 'rho', NF90_DOUBLE, iddim_nrho, id_rho) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_def_var(fgid, 'theta', NF90_DOUBLE, iddim_ntheta, id_theta) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_def_var(fgid, 'xpol', NF90_DOUBLE, (/iddim_nrho, iddim_ntheta /), id_xpol) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_def_var(fgid, 'ypol', NF90_DOUBLE, (/iddim_nrho, iddim_ntheta /), id_ypol) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) ! end of definition nf90_stat = nf90_enddef(fgid) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) ! Compute variables do ip = 1, self%nrho rho_arr(ip) = self%get_rho(ip) enddo do jp = 1, self%ntheta theta_arr(jp) = self%get_theta(jp) enddo do jp = 1, self%ntheta do ip = 1, self%nrho call polar_to_cart(equi, rho_arr(ip), theta_arr(jp), self%phi, xpol(ip, jp), ypol(ip, jp) ) enddo enddo ! Put variables --------------------------------------------------------------------------- nf90_stat = nf90_put_var(fgid, id_rho, rho_arr ) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_put_var(fgid, id_theta, theta_arr ) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_put_var(fgid, id_xpol, xpol ) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_put_var(fgid, id_ypol, ypol ) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) end subroutine module subroutine read_netcdf_polar(self, fgid) class(polar_grid_t), intent(inout) :: self integer, intent(in) :: fgid integer :: nf90_stat, iddim_nrho, iddim_ntheta character(len=*), parameter :: cname = "read_netcdf_polar" ! Name used for logging errors ! Read global attributes (metadata) ------------------------------------------------------- nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'phi', self%phi) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'rhopol_min', self%rhopol_min) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'rhopol_max', self%rhopol_max) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'drho', self%drho) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'dtheta', self%dtheta) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) ! Get dimensions -------------------------------------------------------------------------- nf90_stat = nf90_inq_dimid(fgid, 'nrho', iddim_nrho) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_inquire_dimension(fgid, iddim_nrho, len = self%nrho) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_inq_dimid(fgid, 'ntheta', iddim_ntheta) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_inquire_dimension(fgid, iddim_ntheta, len = self%ntheta) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) end subroutine module subroutine display_polar_grid(self) class(polar_grid_t), intent(in) :: self if (.not. is_master()) then return endif write(get_stdout(), *) & 'Information on polar_grid -----------------------------' write(get_stdout(), 101) & self%phi, & self%rhopol_min, & self%rhopol_max, & self%nrho, & self%ntheta, & self%drho, & self%dtheta 101 FORMAT( 3X,'phi = ',ES14.6E3 /, & 3X,'rhopol_min = ',ES14.6E3 /, & 3X,'rhopol_max = ',ES14.6E3 /, & 3X,'nrho = ',I9 /, & 3X,'ntheta = ',I9 /, & 3X,'drho = ',ES14.6E3 /, & 3X,'dtheta = ',ES14.6E3 ) write(get_stdout(), *) & '-------------------------------------------------------' end subroutine module subroutine destructor(self) type(polar_grid_t), intent(inout) :: self self%phi = FP_NAN self%rhopol_min = FP_NAN self%rhopol_max = FP_NAN self%nrho = -1 self%ntheta = -1 self%drho = FP_NAN self%dtheta = FP_NAN end subroutine end submodule