diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index c0ae90b5ce4d6e0006ad8fd0323754781fe2ded2..527c97121978d31219d156ab88b412ba70bec164 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -53,7 +53,7 @@ MODULE basic
   INTEGER :: lu_stop = 91              ! stop file, see subroutine TESEND
 
   ! To measure computation time
-  real     :: start, finish
+  real(dp) :: start, finish
   real(dp) :: t0_rhs, t0_adv_field, t0_poisson, t0_Sapj, t0_diag, t0_checkfield,&
               t0_step, t0_clos, t0_ghost, t0_coll, t0_process
   real(dp) :: t1_rhs, t1_adv_field, t1_poisson, t1_Sapj, t1_diag, t1_checkfield,&
@@ -165,7 +165,7 @@ CONTAINS
   END SUBROUTINE daytim
   !================================================================================
   SUBROUTINE display_h_min_s(time)
-    real :: time
+    real(dp) :: time
     integer :: days, hours, mins, secs
     days = FLOOR(time/24./3600.);
     hours= FLOOR(time/3600.);
diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
index 72ff5ffa4c9c04a1a1c5d53d7f261b52a2eee15e..9c1f0060ab83aff223aae2935cee35b6bdd9f4ef 100644
--- a/src/calculus_mod.F90
+++ b/src/calculus_mod.F90
@@ -194,7 +194,6 @@ SUBROUTINE simpson_rule_z(f,intf)
    ENDIF
    intf = onethird*deltaz*intf
    ENDIF
-
 END SUBROUTINE simpson_rule_z
 
 END MODULE calculus
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 576e8ae1575819401fb4a35fc344e7ff19caebdd..1160c567169f23462a72d3a2a63f12598d5b05db 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -332,7 +332,7 @@ CONTAINS
 
  SUBROUTINE set_ikx_zBC_map
  IMPLICIT NONE
- REAL :: shift
+ REAL(dp) :: shift
 
  ALLOCATE(ikx_zBC_L(ikys:ikye,ikxs:ikxe))
  ALLOCATE(ikx_zBC_R(ikys:ikye,ikxs:ikxe))
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index cbd94555eb69b43fd261a41b858c9b96470c7606..f2feb2b01d5b828924fb9bcce7a9b559a7f803ed 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -41,7 +41,7 @@ MODULE grid
   REAL(dp), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: zarray
   REAL(dp), DIMENSION(:),   ALLOCATABLE, PUBLIC :: zarray_full
   ! local z weights for computing simpson rule
-  INTEGER,  DIMENSION(:),   ALLOCATABLE, PUBLIC :: zweights_SR
+  REAL(dp),  DIMENSION(:),   ALLOCATABLE, PUBLIC :: zweights_SR
   REAL(dp), PUBLIC, PROTECTED  ::  deltax,  deltay, deltaz, inv_deltaz
   REAL(dp), PUBLIC, PROTECTED  ::  diff_pe_coeff, diff_je_coeff
   REAL(dp), PUBLIC, PROTECTED  ::  diff_pi_coeff, diff_ji_coeff
@@ -373,7 +373,7 @@ CONTAINS
     USE model, ONLY: LINEARITY, N_HD
     IMPLICIT NONE
     REAL(dp), INTENT(IN) :: shear
-    REAL    :: Lx_adapted
+    REAL(dp):: Lx_adapted
     IF(shear .GT. 0) THEN
       IF(my_id.EQ.0) write(*,*) 'Magnetic shear detected: set up sheared kx grid..'
       ! mininal size of box in x to respect dkx = 2pi shear dky
@@ -412,7 +412,7 @@ CONTAINS
         ! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
         local_kxmax = 0._dp
         DO ikx = ikxs,ikxe
-          kxarray(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1)/real(Nkx)))
+          kxarray(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*REAL(ikx-1,dp)/REAL(Nkx,dp)))
           if (ikx .EQ. Nx/2+1)     kxarray(ikx) = -kxarray(ikx)
           ! Finding kx=0
           IF (kxarray(ikx) .EQ. 0) THEN
@@ -429,7 +429,7 @@ CONTAINS
         ! Build the full grids on process 0 to diagnose it without comm
         ! kx
         DO ikx = 1,Nkx
-            kxarray_full(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1)/real(Nkx)))
+            kxarray_full(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*REAL(ikx-1,dp)/REAL(Nkx,dp)))
             IF (ikx .EQ. Nx/2+1) kxarray_full(ikx) = -kxarray_full(ikx)
         END DO
       ELSE ! Odd number of kx (-2 -1 0 1 2)
@@ -493,11 +493,11 @@ CONTAINS
     USE prec_const
     USE model, ONLY: mu_z
     IMPLICIT NONE
-    REAL    :: grid_shift, Lz, zmax, zmin
-    INTEGER :: istart, iend, in
+    REAL(dp) :: grid_shift, Lz, zmax, zmin
+    INTEGER  :: istart, iend, in
     total_nz = Nz
     ! Length of the flux tube (in ballooning angle)
-    Lz         = 2_dp*pi*Npol
+    Lz         = 2_dp*pi*REAL(Npol,dp)
     ! Z stepping (#interval = #points since periodic)
     deltaz        = Lz/REAL(Nz,dp)
     inv_deltaz    = 1._dp/deltaz
@@ -572,11 +572,15 @@ CONTAINS
     ALLOCATE(zweights_SR(izs:ize))
     DO iz = izs,ize
       IF(MODULO(iz,2) .EQ. 1) THEN ! odd iz
-        zweights_SR(iz) = 2._dp
-      ELSE ! even iz
         zweights_SR(iz) = 4._dp
+      ELSE ! even iz
+        zweights_SR(iz) = 2._dp
       ENDIF
     ENDDO
+    ! IF (contains_zmin) &
+    !   zweights_SR(izs) = 1._dp
+    ! IF (contains_zmax) &
+    !   zweights_SR(ize) = 1._dp
   END SUBROUTINE set_zgrid
 
   SUBROUTINE grid_outputinputs(fidres, str)
diff --git a/src/inital.F90 b/src/inital.F90
index f72065887af5a246bf9a94407f0c8e33165a894f..2f61ff9c47e66244605cdff573174a5a8721ad9d 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -306,7 +306,7 @@ SUBROUTINE init_phi
       DO ikx=2,Nkx/2
         phi(iky_0,ikx,izs:ize) = phi(iky_0,Nkx+2-ikx,izs:ize)
       END DO
-      phi(iky_0,ikx_0,izs:ize) = REAL(phi(iky_0,ikx_0,izs:ize)) !origin must be real
+      phi(iky_0,ikx_0,izs:ize) = REAL(phi(iky_0,ikx_0,izs:ize),dp) !origin must be real
     ENDIF
 
     !**** ensure no previous moments initialization
@@ -367,7 +367,7 @@ SUBROUTINE init_phi_ppj
       DO ikx=2,Nkx/2
         phi(iky_0,ikx,izs:ize) = phi(iky_0,Nkx+2-ikx,izs:ize)
       END DO
-      phi(iky_0,ikx_0,izs:ize) = REAL(phi(iky_0,ikx_0,izs:ize)) !origin must be real
+      phi(iky_0,ikx_0,izs:ize) = REAL(phi(iky_0,ikx_0,izs:ize),dp) !origin must be real
     ENDIF
 
     !**** ensure no previous moments initialization
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index e5ccd88c92a06c0a9b70e07e93d4b01f79fbac76..1577957fac953336a3cf6b0819d35a67f92373de 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -98,8 +98,8 @@ SUBROUTINE compute_radial_ion_transport
   pflux_local = integral*iInt_Jacobian
 
   !!!!---------- Sum over all processes ----------
-  buffer(1) = 2._dp*REAL(gflux_local)
-  buffer(2) = 2._dp*REAL(pflux_local)
+  buffer(1) = 2._dp*REAL(gflux_local,dp)
+  buffer(2) = 2._dp*REAL(pflux_local,dp)
   root = 0
   !Gather manually among the rank_p=0 processes and perform the sum
   gflux_ri = 0
@@ -206,8 +206,8 @@ SUBROUTINE compute_radial_electron_transport
   pflux_local = integral*iInt_Jacobian
 
   !!!!---------- Sum over all processes ----------
-  buffer(1) = 2._dp*REAL(gflux_local)
-  buffer(2) = 2._dp*REAL(pflux_local)
+  buffer(1) = 2._dp*REAL(gflux_local,dp)
+  buffer(2) = 2._dp*REAL(pflux_local,dp)
   root = 0
   !Gather manually among the rank_p=0 processes and perform the sum
   gflux_re = 0
@@ -292,7 +292,7 @@ SUBROUTINE compute_radial_ion_heatflux
   call simpson_rule_z(integrant,integral)
   hflux_local = hflux_local + integral*iInt_Jacobian
   ! Double it for taking into account the other half plane
-  buffer(2) = 2._dp*REAL(hflux_local)
+  buffer(2) = 2._dp*REAL(hflux_local,dp)
   root = 0
   !Gather manually among the rank_p=0 processes and perform the sum
   hflux_xi = 0
@@ -375,7 +375,7 @@ SUBROUTINE compute_radial_electron_heatflux
   call simpson_rule_z(integrant,integral)
   hflux_local = hflux_local + integral*iInt_Jacobian
   ! Double it for taking into account the other half plane
-  buffer(2) = 2._dp*REAL(hflux_local)
+  buffer(2) = 2._dp*REAL(hflux_local,dp)
   root = 0
   !Gather manually among the rank_p=0 processes and perform the sum
   hflux_xe = 0
@@ -552,8 +552,8 @@ SUBROUTINE compute_Napjz_spectrum
     DO ikx = ikxs,ikxe
       DO iky = ikys,ikye
         local_sum_e(ips_e:ipe_e,ijs_e:ije_e,izs:ize)  = local_sum_e(ips_e:ipe_e,ijs_e:ije_e,izs:ize)  + &
-        REAL(REAL(moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel)&
-         * CONJG(moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel))))
+        REAL(moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel)&
+         * CONJG(moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel)),dp)
       ENDDO
     ENDDO
     ! sum reduction
diff --git a/src/stepon.F90 b/src/stepon.F90
index b2fec7757aca48fee10ac18350a9d8c56655c1ec..d85f3fd134152d97dacb1a0438df5313833198e3 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -6,6 +6,7 @@ SUBROUTINE stepon
   USE ghosts,                ONLY: update_ghosts_moments, update_ghosts_EM
   use mpi,                   ONLY: MPI_COMM_WORLD
   USE time_integration,      ONLY: ntimelevel
+  USE prec_const,            ONLY: dp
   IMPLICIT NONE
 
   INTEGER :: num_step
@@ -153,7 +154,7 @@ SUBROUTINE stepon
                     moments_e( ip,ij,iky_0,ikx,iz, :) = CONJG(moments_e( ip,ij,iky_0,Nkx+2-ikx,iz, :))
                   END DO
                 ! must be real at origin
-                moments_e(ip,ij, iky_0,ikx_0,iz, :) = REAL(moments_e(ip,ij, iky_0,ikx_0,iz, :))
+                moments_e(ip,ij, iky_0,ikx_0,iz, :) = REAL(moments_e(ip,ij, iky_0,ikx_0,iz, :),dp)
               END DO
             END DO
           END DO
@@ -166,7 +167,7 @@ SUBROUTINE stepon
                   moments_i( ip,ij,iky_0,ikx,iz, :) = CONJG(moments_i( ip,ij,iky_0,Nkx+2-ikx,iz, :))
                 END DO
                 ! must be real at origin and top right
-                moments_i(ip,ij, iky_0,ikx_0,iz, :) = REAL(moments_i(ip,ij, iky_0,ikx_0,iz, :))
+                moments_i(ip,ij, iky_0,ikx_0,iz, :) = REAL(moments_i(ip,ij, iky_0,ikx_0,iz, :),dp)
               END DO
             END DO
           END DO
@@ -175,13 +176,13 @@ SUBROUTINE stepon
             phi(iky_0,ikx,izgs:izge) = phi(iky_0,Nkx+2-ikx,izgs:izge)
           END DO
           ! must be real at origin
-          phi(iky_0,ikx_0,izgs:izge) = REAL(phi(iky_0,ikx_0,izgs:izge))
+          phi(iky_0,ikx_0,izgs:izge) = REAL(phi(iky_0,ikx_0,izgs:izge),dp)
           ! Psi
           DO ikx=2,Nkx/2 !symmetry at ky = 0
             psi(iky_0,ikx,izgs:izge) = psi(iky_0,Nkx+2-ikx,izgs:izge)
           END DO
           ! must be real at origin
-          psi(iky_0,ikx_0,izgs:izge) = REAL(psi(iky_0,ikx_0,izgs:izge))
+          psi(iky_0,ikx_0,izgs:izge) = REAL(psi(iky_0,ikx_0,izgs:izge),dp)
         ENDIF
       END SUBROUTINE enforce_symmetry
 
diff --git a/src/tesend.F90 b/src/tesend.F90
index e908f7e72010f007d405446fab104c6744375422..9a77ea1ceffe90e24445a1acf0298f4f622571c1 100644
--- a/src/tesend.F90
+++ b/src/tesend.F90
@@ -6,7 +6,7 @@ SUBROUTINE tesend
   use prec_const
   IMPLICIT NONE
   LOGICAL :: mlend, mlexist
-  REAL    :: tnow
+  REAL(dp):: tnow
   INTEGER :: ncheck_stop = 100
   CHARACTER(len=*), PARAMETER :: stop_file = 'mystop'