diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 58b4e53c16427cf0638b2770f95a0cfb3e5a8867..72e13ac7064d02aa649f5e550e25bb9834fef4c2 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -77,67 +77,67 @@ contains
   !
   !--------------------------------------------------------------------------------
   ! UNUSED
-  subroutine eval_salpha_geometry
-   ! evaluate s-alpha geometry model
-   implicit none
-   REAL(dp) :: z, kx, ky
-
-   zloop: DO iz = izs,ize
-    z = zarray(iz)
-    gxx(iz) = 1._dp
-    gxy(iz) = shear*z
-    gyy(iz) = 1._dp + (shear*z)**2
-    gyz(iz) = 1._dp/eps
-    gxz(iz) = 0._dp
-
-   ! Relative strengh of radius
-      hatR(iz) = 1._dp + eps*cos(z)
-
-   ! Jacobian
-      Jacobian(iz) = q0*hatR(iz)
-
-   ! Relative strengh of modulus of B
-      hatB(iz) = 1._dp / hatR( iz)
-
-   ! Derivative of the magnetic field strenght
-      gradxB(iz) = - cos( z) / hatR(iz)
-      gradzB(iz) = eps * sin(z)
-
-   ! Gemoetrical coefficients for the curvature operator
-   ! Note: Gamma2 and Gamma3 are obtained directly form Gamma1 in the expression of the curvature operator implemented here
-   !
-      Gamma1(iz) = gxy(iz) * gxy(iz) - gxx(iz) * gyy(iz)
-      Gamma2(iz) = gxz(iz) * gxy(iz) - gxx(iz) * gyz(iz)
-      Gamma3(iz) = gxz(iz) * gyy(iz) - gxy(iz) * gyz(iz)
-
-   ! Curvature operator
-      DO iky = ikys, ikye
-        ky = kyarray(iky)
-         DO ikx= ikxl, ikxr
-           kx = kxarray(ikx,iky)
-           Cxy(ikx, iky, iz) = (-sin(z)*kx - (cos(z) + shear* z* sin(z))*ky) * hatB(iz) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
-         ENDDO
-      ENDDO
-
-   ! coefficient in the front of parallel derivative
-      gradz_coeff(iz) = 1._dp / Jacobian(iz) / hatB(iz)
-   ENDDO
-
-   ! Evaluate perpendicular wavenumber
-   !  k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy}
-   !  normalized to rhos_
-      DO iky = ikys, ikye
-        ky = kyarray(iky)
-         DO ikx = ikxl, ikxr
-           kx = kxarray(ikx,iky)
-            kperp_array(ikx, iky, iz) = sqrt( gxx(iz)*kx**2 + 2._dp* gxy(iz) * kx*ky + gyy(iz)* ky**2) !! / hatB( iz) ! there is a factor 1/B from the normalization; important to match GENE
-         ENDDO
-      ENDDO
-   ENDDO zloop
-   !
-   IF( me .eq. 0 ) PRINT*, 'max kperp = ', maxval( kperp_array)
-
-  END SUBROUTINE eval_salpha_geometry
+  ! subroutine eval_salpha_geometry
+  !  ! evaluate s-alpha geometry model
+  !  implicit none
+  !  REAL(dp) :: z, kx, ky
+  !
+  !  zloop: DO iz = izs,ize
+  !   z = zarray(iz)
+  !   gxx(iz) = 1._dp
+  !   gxy(iz) = shear*z
+  !   gyy(iz) = 1._dp + (shear*z)**2
+  !   gyz(iz) = 1._dp/eps
+  !   gxz(iz) = 0._dp
+  !
+  !  ! Relative strengh of radius
+  !     hatR(iz) = 1._dp + eps*cos(z)
+  !
+  !  ! Jacobian
+  !     Jacobian(iz) = q0*hatR(iz)
+  !
+  !  ! Relative strengh of modulus of B
+  !     hatB(iz) = 1._dp / hatR( iz)
+  !
+  !  ! Derivative of the magnetic field strenght
+  !     gradxB(iz) = - cos( z) / hatR(iz)
+  !     gradzB(iz) = eps * sin(z)
+  !
+  !  ! Gemoetrical coefficients for the curvature operator
+  !  ! Note: Gamma2 and Gamma3 are obtained directly form Gamma1 in the expression of the curvature operator implemented here
+  !  !
+  !     Gamma1(iz) = gxy(iz) * gxy(iz) - gxx(iz) * gyy(iz)
+  !     Gamma2(iz) = gxz(iz) * gxy(iz) - gxx(iz) * gyz(iz)
+  !     Gamma3(iz) = gxz(iz) * gyy(iz) - gxy(iz) * gyz(iz)
+  !
+  !  ! Curvature operator
+  !     DO iky = ikys, ikye
+  !       ky = kyarray(iky)
+  !        DO ikx= ikxl, ikxr
+  !          kx = kxarray(ikx,iky)
+  !          Cxy(ikx, iky, iz) = (-sin(z)*kx - (cos(z) + shear* z* sin(z))*ky) * hatB(iz) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+  !        ENDDO
+  !     ENDDO
+  !
+  !  ! coefficient in the front of parallel derivative
+  !     gradz_coeff(iz) = 1._dp / Jacobian(iz) / hatB(iz)
+  !  ENDDO
+  !
+  !  ! Evaluate perpendicular wavenumber
+  !  !  k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy}
+  !  !  normalized to rhos_
+  !     DO iky = ikys, ikye
+  !       ky = kyarray(iky)
+  !        DO ikx = ikxl, ikxr
+  !          kx = kxarray(ikx,iky)
+  !           kperp_array(ikx, iky, iz) = sqrt( gxx(iz)*kx**2 + 2._dp* gxy(iz) * kx*ky + gyy(iz)* ky**2) !! / hatB( iz) ! there is a factor 1/B from the normalization; important to match GENE
+  !        ENDDO
+  !     ENDDO
+  !  ENDDO zloop
+  !  !
+  !  IF( me .eq. 0 ) PRINT*, 'max kperp = ', maxval( kperp_array)
+  !
+  ! END SUBROUTINE eval_salpha_geometry
   !
   !--------------------------------------------------------------------------------
   !