`btor` Initialization Missing for VMEC Geometries, Causing NaNs in Momentum Flux Calculations

Issue Description

Currently, btor, used for calculating the momentum flux, is not initialized or computed for VMEC geometries. This leads to NaNs and incorrect results in the output.

Code Example

See stella_geometry.f90 line 300:

case (geo_option_vmec)
    vmec_chosen = .true.
    !> read in input parameters for vmec
    !> nalpha may be specified via input file
    if (debug) write (*, *) 'init_geometry::read_vmec_parameters'
    call read_vmec_parameters
    !> allocate geometry arrays
    if (debug) write (*, *) 'init_geometry::allocate_arrays'
    call allocate_arrays(nalpha, nzgrid)
    if (debug) write (*, *) 'init_geometry::allocate_temporary_arrays'
    call allocate_temporary_arrays(nalpha, nzgrid)
    !> get geometry coefficients from vmec
    if (debug) write (*, *) 'init_geometry::get_vmec_geo'
    !> abs(twist_and_shift_geo_fac) is dkx/dky * jtwist
    !> minus its sign gives the direction of the shift in kx
    !> to be used for twist-and-shift BC
    allocate (twist_and_shift_geo_fac_full(nalpha, -nzgrid:nzgrid))
    twist_and_shift_geo_fac_full = 0
    !> The code will start computing for the initial zgrid
    stellarator_symmetric_BC = .false.
    call get_vmec_geo(new_zeta_min, stellarator_symmetric_BC, &
                      nzgrid, nalpha, naky, geo_surf, grho, bmag, gradpar, &
                      b_dot_grad_z, grad_alpha_grad_alpha, &
                      grad_alpha_grad_psi, grad_psi_grad_psi, &
                      gds23, gds24, gds25, gds26, gbdrift_alpha, gbdrift0_psi, &
                      cvdrift_alpha, cvdrift0_psi, sign_torflux, &
                      theta_vmec, zed_scalefac, aref, bref, alpha, zeta, &
                      field_period_ratio, x_displacement_fac)

See stella_diagnostics.f90 line 577:

! get momentum flux
! parallel component
gtmp1 = g(:, :, ikxkyz) * spread(vpa, 2, nmu) * geo_surf%rmaj * btor(iz) / bmag(ia, iz)
call gyro_average(gtmp1, ikxkyz, gtmp2)
gtmp1 = -g(:, :, ikxkyz) * zi * aky(iky) * spread(vperp2(ia, iz, :), 1, nvpa) * geo_surf%rhoc &
        * (gds21(ia, iz) + theta0(iky, ikx) * gds22(ia, iz)) * spec(is)%smz &
        / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2)
call gyro_average_j1(gtmp1, ikxkyz, gtmp3)
gtmp1 = gtmp2 + gtmp3

call get_one_flux(iky, iz, flx_norm(iz), gtmp1, phi(iky, ikx, iz, it), vflx(is))
call get_one_flux(iky, iz, flx_norm_partial, gtmp1, phi(iky, ikx, iz, it), vflx_vs_kxkyz(iky, ikx, iz, it, is))

Expected Behavior

btor should be properly initialized and computed for VMEC geometries to avoid NaNs and ensure correct results.

stellaGK / stella

`btor` Initialization Missing for VMEC Geometries, Causing NaNs in Momentum Flux Calculations #131