Fix warnings and stabilize orbit diagnostics

app/simple_diag_orbit.f90

@@ -19,15 +19,11 @@ program diag_orbit_main
 
 implicit none
 
-! Define dp kind parameter
-integer, parameter :: dp = kind(1.0d0)
-
 character(256) :: config_file, particle_arg
 type(tracer_t) :: norb
 type(symplectic_integrator_t) :: si
 type(field_can_t) :: field_can
 integer :: particle_number
-real(dp), dimension(5) :: z
 
 ! Initialize timing
 call init_timer()
@@ -143,4 +139,4 @@ program diag_orbit_main
 print *, "Generated detailed trajectory plots showing real particle"
 print *, "motion in the magnetic field using symplectic integration."
 
-end program diag_orbit_main
+end program diag_orbit_main

cmake/Util.cmake

@@ -61,14 +61,6 @@ function(find_or_fetch DEPENDENCY)
         endif()
     endif()
 
-    if(NOT _use_local AND DEFINED ENV{CODE})
-        get_filename_component(_candidate "$ENV{CODE}/${_dep_lower}" ABSOLUTE)
-        if(IS_DIRECTORY "${_candidate}")
-            set(_source_dir "${_candidate}")
-            set(_use_local TRUE)
-        endif()
-    endif()
-
     if(_use_local)
         if("${_dep_lower}" STREQUAL "libneo")
             set(LIBNEO_ENABLE_TESTS OFF CACHE BOOL "" FORCE)

python/CMakeLists.txt

@@ -86,11 +86,7 @@ set(FILES_TO_WRAP
 )
 
 # Add libneo files to wrap
-if(DEFINED ENV{CODE})
-    set(LIBNEO_SOURCE_DIR $ENV{CODE}/libneo)
-else()
-    set(LIBNEO_SOURCE_DIR ${libneo_SOURCE_DIR})
-endif()
+set(LIBNEO_SOURCE_DIR ${libneo_SOURCE_DIR})
 
 list(APPEND FILES_TO_WRAP ${LIBNEO_SOURCE_DIR}/src/canonical_coordinates_mod.f90)
 

src/boozer_converter.F90

@@ -82,7 +82,7 @@ subroutine get_boozer_coordinates_impl
         use vector_potentail_mod, only: ns, hs
         use new_vmec_stuff_mod, only: n_theta, n_phi, h_theta, h_phi, ns_s, ns_tp
         use boozer_coordinates_mod, only: ns_s_B, ns_tp_B, ns_B, n_theta_B, n_phi_B, &
-                                          hs_B, h_theta_B, h_phi_B, use_B_r
+                                          hs_B, h_theta_B, h_phi_B
 
         implicit none
 
@@ -784,7 +784,7 @@ end subroutine compute_boozer_data
 
     !> Compute radial covariant magnetic field B_rho from symmetry flux coordinates
     subroutine compute_br_from_symflux(rho_tor, aiota_arr, Gfunc, Bcovar_symfl)
-        use boozer_coordinates_mod, only: ns_B, n_theta_B, n_phi_B
+        use boozer_coordinates_mod, only: ns_B, n_phi_B
         use plag_coeff_sub, only: plag_coeff
 
         real(dp), intent(in) :: rho_tor(:)
@@ -798,7 +798,7 @@ subroutine compute_br_from_symflux(rho_tor, aiota_arr, Gfunc, Bcovar_symfl)
         integer :: i_rho, i_phi, ibeg, iend, nshift
         real(dp) :: coef(0:NDER, NPOILAG)
 
-        nshift = NPOILAG/2
+        nshift = (NPOILAG - 1)/2
 
         do i_rho = 1, ns_B
             ibeg = i_rho - nshift

src/classification.f90

@@ -1,4 +1,5 @@
 module classification
+use, intrinsic :: iso_fortran_env, only: dp => real64
 use omp_lib
 use params, only: zstart, zend, times_lost, trap_par, perp_inv, iclass, &
     ntimstep, confpart_trap, confpart_pass, notrace_passing, contr_pp, &
@@ -18,9 +19,6 @@ module classification
 
   implicit none
 
-  ! Define real(dp) kind parameter
-  integer, parameter :: dp = kind(1.0d0)
-
   ! Classification result type - separates data from I/O
   ! Note: 0=unclassified means the classification was not computed
   ! This depends on orbit type (trapped/passing) and class_plot flag
@@ -65,7 +63,7 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result)
     type(tracer_t), intent(inout) :: anorb
     integer, intent(in) :: ipart
     type(classification_result_t), intent(out) :: class_result
-    integer :: ierr, ierr_coll
+    integer :: ierr
     real(dp), dimension(5) :: z
     real(dp) :: bmod,sqrtg
     real(dp), dimension(3) :: bder, hcovar, hctrvr, hcurl

src/cut_detector.f90

@@ -1,4 +1,5 @@
 module cut_detector
+    use, intrinsic :: iso_fortran_env, only: dp => real64
     use util, only: twopi
     use simple, only: tstep
     use orbit_symplectic, only: symplectic_integrator_t
@@ -9,8 +10,6 @@ module cut_detector
 
     implicit none
 
-    ! Define real(dp) kind parameter
-    integer, parameter :: dp = kind(1.0d0)
     save
 
     integer, parameter :: n_tip_vars = 6
@@ -86,33 +85,29 @@ subroutine trace_to_cut(self, si, f, z, var_cut, cut_type, ierr)
       integer :: i
       real(dp) :: phiper = 0.0d0
 
-      do i=1, nstep_max
+      do i = 1, nplagr
         call tstep(si, f, z, ierr)
-        if(ierr.ne.0) exit
-
-        self%par_inv = self%par_inv+z(5)**2 ! parallel adiabatic invariant
-
-        if(i.le.nplagr) then          !<=first nplagr points to initialize stencil
-          ! Note: i is guaranteed to be <= nplagr here, so array access is safe
-          self%orb_sten(1:5,i)=z
-          self%orb_sten(6,i)=self%par_inv
-        else                          !<=normal case, shift stencil
-          self%orb_sten(1:5,self%ipoi(1))=z
-          self%orb_sten(6,self%ipoi(1))=self%par_inv
-          self%ipoi=cshift(self%ipoi,1)
-        endif
+        if (ierr.ne.0) exit
+
+        self%par_inv = self%par_inv + z(5)**2 ! parallel adiabatic invariant
+
+        self%orb_sten(1:5, i) = z
+        self%orb_sten(6, i) = self%par_inv
 
         !-------------------------------------------------------------------------
         ! Tip detection and interpolation
 
-        if(self%alam_prev.lt.0.d0.and.z(5).gt.0.d0) self%itip=0   !<=tip has been passed
+        if(self%alam_prev.lt.0.d0.and.z(5).gt.0.d0) then
+          self%itip=0   !<=tip has been passed
+        end if
         self%itip=self%itip+1
         self%alam_prev=z(5)
 
         !<=use only initialized stencil
         if(i.gt.nplagr .and. self%itip.eq.nplagr/2) then
           cut_type = 0
-          call plag_coeff(nplagr, nder, 0d0, self%orb_sten(5, self%ipoi), self%coef)
+          call plag_coeff(nplagr, nder, 0d0, self%orb_sten(5, self%ipoi), &
+            self%coef)
           var_cut = matmul(self%orb_sten(:, self%ipoi), self%coef(0,:))
           var_cut(2) = modulo(var_cut(2), twopi)
           var_cut(3) = modulo(var_cut(3), twopi)
@@ -140,7 +135,70 @@ subroutine trace_to_cut(self, si, f, z, var_cut, cut_type, ierr)
         if(i.gt.nplagr .and. self%iper.eq.nplagr/2) then
           cut_type = 1
           !<=stencil around periodic boundary is complete, interpolate
-          call plag_coeff(nplagr, nder, 0d0, self%orb_sten(3, self%ipoi) - phiper, self%coef)
+          call plag_coeff(nplagr, nder, 0d0, self%orb_sten(3, self%ipoi) - &
+            phiper, self%coef)
+          var_cut = matmul(self%orb_sten(:, self%ipoi), self%coef(0,:))
+          var_cut(2) = modulo(var_cut(2), twopi)
+          var_cut(3) = modulo(var_cut(3), twopi)
+          return
+        end if
+
+        ! End periodic boundary footprint detection and interpolation
+        !-------------------------------------------------------------------------
+
+      end do
+
+      do i = nplagr + 1, nstep_max
+        call tstep(si, f, z, ierr)
+        if (ierr.ne.0) exit
+
+        self%par_inv = self%par_inv + z(5)**2 ! parallel adiabatic invariant
+
+        self%orb_sten(1:5, self%ipoi(1)) = z
+        self%orb_sten(6, self%ipoi(1)) = self%par_inv
+        self%ipoi = cshift(self%ipoi, 1)
+
+        !-------------------------------------------------------------------------
+        ! Tip detection and interpolation
+
+        if(self%alam_prev.lt.0.d0.and.z(5).gt.0.d0) then
+          self%itip=0   !<=tip has been passed
+        end if
+        self%itip=self%itip+1
+        self%alam_prev=z(5)
+
+        if(self%itip.eq.nplagr/2) then
+          cut_type = 0
+          call plag_coeff(nplagr, nder, 0d0, self%orb_sten(5, self%ipoi), &
+            self%coef)
+          var_cut = matmul(self%orb_sten(:, self%ipoi), self%coef(0,:))
+          var_cut(2) = modulo(var_cut(2), twopi)
+          var_cut(3) = modulo(var_cut(3), twopi)
+          self%par_inv = self%par_inv - var_cut(6)
+          return
+        end if
+
+        ! End tip detection and interpolation
+
+        !-------------------------------------------------------------------------
+        ! Periodic boundary footprint detection and interpolation
+
+        if(z(3).gt.dble(self%kper+1)*self%fper) then
+          self%iper=0   !<=periodic boundary has been passed
+          phiper=dble(self%kper+1)*self%fper
+          self%kper=self%kper+1
+        elseif(z(3).lt.dble(self%kper)*self%fper) then
+          self%iper=0   !<=periodic boundary has been passed
+          phiper=dble(self%kper)*self%fper
+          self%kper=self%kper-1
+        endif
+        self%iper=self%iper+1
+
+        if(self%iper.eq.nplagr/2) then
+          cut_type = 1
+          !<=stencil around periodic boundary is complete, interpolate
+          call plag_coeff(nplagr, nder, 0d0, self%orb_sten(3, self%ipoi) - &
+            phiper, self%coef)
           var_cut = matmul(self%orb_sten(:, self%ipoi), self%coef(0,:))
           var_cut(2) = modulo(var_cut(2), twopi)
           var_cut(3) = modulo(var_cut(3), twopi)

src/diag/diag_orbit.f90

@@ -133,7 +133,7 @@ subroutine integrate_orbit_with_trajectory_debug(si, f, particle_number)
     real(dp), allocatable :: pphi_traj(:)
     integer, allocatable :: newton_iter_traj(:)
     real(dp), dimension(5) :: z, xlast
-    integer :: it, ktau, point_idx, newton_iters, ierr_orbit
+    integer :: it, ktau, point_idx, newton_iters
     integer, parameter :: maxit = 32
     real(dp) :: current_time
     integer :: total_points, field_eval_count
@@ -269,4 +269,4 @@ subroutine write_trajectory_data(time_traj, s_traj, theta_traj, phi_traj, pphi_t
 
 end subroutine write_trajectory_data
 
-end module diag_orbit
+end module diag_orbit

src/netcdf_results_output.f90

@@ -318,7 +318,7 @@ subroutine compute_cartesian_positions(xstart_cart, xend_cart)
 
             ! Check if particle was traced: zend is set to exactly 0 for untraced
             ! particles (see simple_main.f90:493 and classification.f90:94)
-            zend_is_zero = all(zend(1:3, i) == 0.0_dp)
+            zend_is_zero = all(abs(zend(1:3, i)) <= 0.0_dp)
             if (zend_is_zero) then
                 ! Particle not traced - use start position
                 xend_cart(:, i) = xstart_cart(:, i)

src/orbit_symplectic.f90

@@ -1,5 +1,6 @@
 module orbit_symplectic
 
+use, intrinsic :: iso_fortran_env, only: dp => real64
 use util, only: pi, twopi
 use field_can_mod, only: field_can_t, get_val, get_derivatives, get_derivatives2, &
   eval_field => evaluate
@@ -15,9 +16,6 @@ module orbit_symplectic
 
 implicit none
 
-! Define real(dp) kind parameter
-integer, parameter :: dp = kind(1.0d0)
-
 procedure(orbit_timestep_sympl_i), pointer :: orbit_timestep_sympl => null()
 
 contains
@@ -26,8 +24,6 @@ module orbit_symplectic
   !
 subroutine orbit_sympl_init(si, f, z, dt, ntau, rtol_init, mode_init)
   !
-  use plag_coeff_sub, only : plag_coeff
-
   type(symplectic_integrator_t), intent(inout) :: si
   type(field_can_t), intent(inout) :: f
   real(dp), intent(in) :: z(:)
@@ -36,8 +32,6 @@ subroutine orbit_sympl_init(si, f, z, dt, ntau, rtol_init, mode_init)
   real(dp), intent(in) :: rtol_init
   integer, intent(in) :: mode_init
 
-  integer :: k
-
   si%atol = 1d-15
   si%rtol = rtol_init
 
@@ -429,7 +423,6 @@ subroutine newton2(si, f, x, atol, rtol, maxit, xlast)
   real(dp), intent(out) :: xlast(n)
 
   real(dp) :: fvec(n), fjac(n,n), jinv(n,n)
-  integer :: pivot(n), info
 
   real(dp) :: xabs(n), tolref(n), fabs(n)
   real(dp) :: det
@@ -792,8 +785,13 @@ subroutine jac_rk_lobatto(si, fs, s, jac)
 
   call coeff_rk_lobatto(s, a, ahat, b, c)
   jac = 0d0
+  dHprime = 0.0d0
 
-  Hprime = fs%dH(1)/fs%dpth(1)
+  Hprime = 0.0d0
+  Hprime(1) = fs(1)%dH(1)/fs(1)%dpth(1)
+  do k = 2, s
+    Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1)
+  end do
   dHprime(1) = (fs(1)%d2H(1)-Hprime(1)*fs(1)%d2pth(1))/fs(1)%dpth(1)  ! d/dr
   dHprime(2) = (fs(1)%d2H(7)-Hprime(1)*fs(1)%d2pth(7))/fs(1)%dpth(1)  ! d/dpph
   do k = 2, s
@@ -1225,7 +1223,7 @@ subroutine orbit_timestep_sympl_expl_impl_euler(si, f, ierr)
   integer, parameter :: maxit = 32
 
   real(dp), dimension(n) :: x, xlast
-  integer :: k, ktau
+  integer :: ktau
 
   ierr = 0
   ktau = 0
@@ -1284,7 +1282,7 @@ subroutine orbit_timestep_sympl_impl_expl_euler(si, f, ierr)
   integer, parameter :: maxit = 32
 
   real(dp), dimension(n) :: x, xlast, dz
-  integer :: k, ktau
+  integer :: ktau
 
   ierr = 0
   ktau = 0
@@ -1346,7 +1344,7 @@ subroutine orbit_timestep_sympl_midpoint(si, f, ierr)
   integer, parameter :: maxit = 8
 
   real(dp), dimension(n) :: x, xlast
-  integer :: k, ktau
+  integer :: ktau
 
   ierr = 0
   ktau = 0

src/samplers.f90

@@ -1,11 +1,9 @@
 module samplers
+  use, intrinsic :: iso_fortran_env, only: dp => real64
   use util
 
   implicit none
 
-  ! Define real(dp) kind parameter
-  integer, parameter :: dp = kind(1.0d0)
-
   character(len=*), parameter :: START_FILE = 'start.dat'
   character(len=*), parameter :: START_FILE_ANTS = 'start_ants.dat'
   character(len=*), parameter :: START_FILE_BATCH = 'batch.dat'
@@ -158,7 +156,7 @@ end subroutine sample_surface_fieldline
 
   subroutine sample_grid(zstart, grid_density)
     use params, only: ntestpart, zstart_dim1, zend, times_lost, &
-        trap_par, perp_inv, iclass, xstart, sbeg
+        trap_par, perp_inv, iclass, sbeg
     use util, only: pi
 
     real(dp), dimension(:,:), allocatable, intent(inout) :: zstart