diff --git a/fort.90 b/fort.90
index 88f8dbfb548370a7d1f198fd10a33190de3db03d..3638ba08a70dd3b4ab33541847a86b5341f6537d 100644
--- a/fort.90
+++ b/fort.90
@@ -19,8 +19,12 @@
 &GEOMETRY
   geom   = 'Z-pinch'
   q0     = 0.0
-  shear  = 0.0 
+  shear  = 0.0
   eps    = 0.0
+  !geom   = 'miller'
+  !q0     = 1.4
+  !shear  = 0.8
+  !eps    = 0.18
   kappa  = 1.0
   s_kappa= 0.0
   delta  = 0.0
@@ -29,6 +33,7 @@
   s_zeta = 0.0
   parallel_bc = 'shearless'
   shift_y= 0.0
+  Npol   = 1.0
 /
 &OUTPUT_PAR
   dtsave_0d = 1
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index ffcedfdc2081b3eaae89ea28c574be1ad9f7dbd2..552a48c97d1764e410789f315936fe6d87fe788c 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -20,7 +20,7 @@ implicit none
   REAL(xp),    PUBLIC, PROTECTED :: s_zeta    = 0._xp ! '' szeta = r dzeta/dr
   ! to apply shift in the parallel z-BC if shearless
   REAL(xp),    PUBLIC, PROTECTED :: shift_y   = 0._xp ! for Arno <3
-  INTEGER,     PUBLIC, PROTECTED :: Npol      = 1         ! number of poloidal turns
+  REAL(xp),    PUBLIC, PROTECTED :: Npol      = 1._xp ! number of poloidal turns (real for 3D zpinch studies)
   ! Chooses the type of parallel BC we use for the unconnected kx modes (active for non-zero shear only)
   !  'periodic'     : Connect a disconnected kx to a mode on the other cadran
   !  'dirichlet'    : Connect a disconnected kx to 0
@@ -103,7 +103,7 @@ CONTAINS
     REAL(xp) :: kx,ky
     COMPLEX(xp), DIMENSION(local_nz) :: integrant
     real(xp) :: G1,G2,G3,Cx,Cy
-    INTEGER  :: eo,iz,iky,ikx
+    INTEGER  :: eo,iz,iky,ikx, Npol_i
 
     ! Allocate arrays
     CALL geometry_allocate_mem(local_nky,local_nkx,local_nz,ngz,nzgrid)
@@ -114,6 +114,7 @@ CONTAINS
     ELSE
       SELECT CASE(geom)
         CASE('s-alpha')
+          IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for s-alpha geometry"
           CALL speak('s-alpha geometry')
           call eval_salpha_geometry
         CASE('Z-pinch','z-pinch','Zpinch','zpinch')
@@ -126,8 +127,9 @@ CONTAINS
           kappa   = 1._xp
         CASE('miller')
           CALL speak('Miller geometry')
+          IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for Miller geometry"
           call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta)
-          call get_miller(eps,major_R,major_Z,q0,shear,Npol,alpha_MHD,edge_opt,&
+          call get_miller(eps,major_R,major_Z,q0,shear,FLOOR(Npol),alpha_MHD,edge_opt,&
                           C_y,C_xy,xpdx_pm_geom,gxx,gyy,gzz,gxy,gxz,gyz,&
                           dBdx,dBdy,hatB,jacobian,dBdz,hatR,hatZ,dxdR,dxdZ,&
                           Ckxky,gradz_coeff)
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 4817c7b84040025ce0aa01558defc117478a6329..89f533248a7f2d6624f7c2263d5ccc5fd65360f1 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -146,8 +146,7 @@ CONTAINS
   !  1 2 3 4
   SUBROUTINE set_grids(shear,Npol,LINEARITY,N_HD,EM,Na)
     USE fourier, ONLY: init_grid_distr_and_plans
-    REAL(xp), INTENT(IN) :: shear
-    INTEGER,  INTENT(IN) :: Npol
+    REAL(xp), INTENT(IN) :: shear, Npol
     CHARACTER(len=*), INTENT(IN) :: LINEARITY
     INTEGER, INTENT(IN)  :: N_HD
     LOGICAL, INTENT(IN)  :: EM
@@ -398,8 +397,7 @@ CONTAINS
   SUBROUTINE set_kxgrid(shear,Npol,LINEARITY,N_HD)
     USE prec_const
     IMPLICIT NONE
-    REAL(xp), INTENT(IN) :: shear
-    INTEGER,  INTENT(IN) :: Npol
+    REAL(xp), INTENT(IN) :: shear, Npol
     CHARACTER(len=*), INTENT(IN) ::LINEARITY
     INTEGER, INTENT(IN)  :: N_HD
     INTEGER :: ikx
@@ -490,11 +488,11 @@ CONTAINS
     USE prec_const
     USE parallel, ONLY: num_procs_z, rank_z
     IMPLICIT NONE
-    REAL(xp):: grid_shift, Lz, zmax, zmin
-    INTEGER :: istart, iend, in, Npol, iz, ig, eo, iglob
+    REAL(xp):: grid_shift, Lz, zmax, zmin, Npol
+    INTEGER :: istart, iend, in, iz, ig, eo, iglob
     total_nz = Nz
     ! Length of the flux tube (in ballooning angle)
-    Lz         = 2._xp*pi*REAL(Npol,xp)
+    Lz         = 2._xp*pi*Npol
     ! Z stepping (#interval = #points since periodic)
     deltaz        = Lz/REAL(Nz,xp)
     inv_deltaz    = 1._xp/deltaz
@@ -515,24 +513,24 @@ CONTAINS
     ENDIF
     ! Build the full grids on process 0 to diagnose it without comm
     ALLOCATE(zarray_full(total_nz))
-    IF (Nz .EQ. 1) Npol = 0
-    zmax = 0; zmin = 0;
-    DO iz = 1,total_nz ! z in [-pi pi-dz] x Npol
-      zarray_full(iz) = REAL(iz-1,xp)*deltaz - Lz/2._xp
-      IF(zarray_full(iz) .GT. zmax) zmax = zarray_full(iz)
-      IF(zarray_full(iz) .LT. zmin) zmin = zarray_full(iz)
-    END DO
-    !! Parallel data distribution
-    IF( (Nz .EQ. 1) .AND. (num_procs_z .GT. 1) ) &
-    ERROR STOP '>> ERROR << Cannot have multiple core in z-direction (Nz = 1)'
-    ! Local data distribution
-    CALL decomp1D(total_nz, num_procs_z, rank_z, izs, ize)
-    local_nz        = ize - izs + 1
-    local_nz_offset = izs - 1
-    ! Ghosts boundaries (depend on the order of z operators)
-    IF(Nz .EQ. 1) THEN
-      ngz              = 0
-      zarray_full(izs) = 0
+    IF (Nz .EQ. 1) Npol = 0._xp
+      zmax = 0; zmin = 0;
+      DO iz = 1,total_nz ! z in [-pi pi-dz] x Npol
+        zarray_full(iz) = REAL(iz-1,xp)*deltaz - Lz/2._xp
+        IF(zarray_full(iz) .GT. zmax) zmax = zarray_full(iz)
+        IF(zarray_full(iz) .LT. zmin) zmin = zarray_full(iz)
+      END DO
+      !! Parallel data distribution
+      IF( (Nz .EQ. 1) .AND. (num_procs_z .GT. 1) ) &
+      ERROR STOP '>> ERROR << Cannot have multiple core in z-direction (Nz = 1)'
+      ! Local data distribution
+      CALL decomp1D(total_nz, num_procs_z, rank_z, izs, ize)
+      local_nz        = ize - izs + 1
+      local_nz_offset = izs - 1
+      ! Ghosts boundaries (depend on the order of z operators)
+      IF(Nz .EQ. 1) THEN
+        ngz              = 0
+        zarray_full(izs) = 0
     ELSEIF(Nz .GE. 4) THEN
       ngz =4
       IF(mod(Nz,2) .NE. 0 ) THEN