diff --git a/src/compute_gradients.F90 b/src/compute_gradients.F90
new file mode 100644
index 0000000000000000000000000000000000000000..06c920d688846c9a0af11ece0e031cf05e4295a3
--- /dev/null
+++ b/src/compute_gradients.F90
@@ -0,0 +1,28 @@
+SUBROUTINE compute_gradients
+    !   This routine compute the gradients of the moments without copying or slices.
+    ! It should be faster than using a routine taking a slice as argument (see 
+    ! calculus_mod) since it avoid copying.
+    USE fields, ONLY: moments
+    USE grid,   ONLY: local_nz, ngz, inv_deltaz
+    USE prec_const, ONLY: dp
+    IMPLICIT NONE
+
+    REAL(dp), dimension(-2:2) :: dz_usu = &
+     (/  1._dp/12._dp, -2._dp/3._dp, 0._dp, 2._dp/3._dp, -1._dp/12._dp /) ! fd4 centered stencil
+    REAL(dp), dimension(-2:1) :: dz_o2e = &
+     (/ 1._dp/24._dp,-9._dp/8._dp, 9._dp/8._dp,-1._dp/24._dp /) ! fd4 odd to even stencil
+    REAL(dp), dimension(-1:2) :: dz_e2o = &
+     (/ 1._dp/24._dp,-9._dp/8._dp, 9._dp/8._dp,-1._dp/24._dp /) ! fd4 odd to even stencil
+     REAL(dp), dimension(-2:2) :: dz2_usu = &
+     (/-1._dp/12._dp, 4._dp/3._dp, -5._dp/2._dp, 4._dp/3._dp, -1._dp/12._dp /)! 2th derivative, 4th order (for parallel hypdiff)
+     REAL(dp), dimension(-2:2) :: dz4_usu = &
+     (/  1._dp, -4._dp, 6._dp, -4._dp, 1._dp /) ! 4th derivative, 2nd order (for parallel hypdiff)
+     REAL(dp), dimension(-2:1) :: iz_o2e = &
+     (/ -1._dp/16._dp, 9._dp/16._dp, 9._dp/16._dp, -1._dp/16._dp /) ! grid interpolation, 4th order, odd to even
+     REAL(dp), dimension(-1:2) :: iz_e2o = &
+     (/ -1._dp/16._dp, 9._dp/16._dp, 9._dp/16._dp, -1._dp/16._dp /) ! grid interpolation, 4th order, even to odd
+    PUBLIC :: simpson_rule_z, interp_z, grad_z, grad_z4
+  
+IMPLICIT NONE
+
+END SUBROUTINE compute_gradients
\ No newline at end of file
diff --git a/testcases/cyclone_example/fort_00.90 b/testcases/cyclone_example/fort_00.90
new file mode 100644
index 0000000000000000000000000000000000000000..44ccd9d7a4c9cb1aeb84bae059fce1ade6ee1757
--- /dev/null
+++ b/testcases/cyclone_example/fort_00.90
@@ -0,0 +1,92 @@
+&BASIC
+  nrun       = 99999999
+  dt         = 0.01
+  tmax       = 50
+  maxruntime = 356400
+  job2load   = -1
+/
+&GRID
+  pmax   = 4
+  jmax   = 2
+  Nx     = 128
+  Lx     = 120
+  Ny     = 48
+  Ly     = 120
+  Nz     = 16
+  SG     = .f.
+  Nexc   = 0
+/
+&GEOMETRY
+  geom   = 's-alpha'
+  q0     = 1.4
+  shear  = 0.8
+  eps    = 0.18
+  kappa  = 1.0
+  s_kappa= 0.0
+  delta  = 0.0
+  s_delta= 0.0
+  zeta   = 0.0
+  s_zeta = 0.0
+  parallel_bc = 'dirichlet'
+  shift_y= 0.0
+/
+&OUTPUT_PAR
+  dtsave_0d = 0.1
+  dtsave_1d = -1
+  dtsave_2d = -1
+  dtsave_3d = 1
+  dtsave_5d = 10
+  write_doubleprecision = .f.
+  write_gamma = .t.
+  write_hf    = .t.
+  write_phi   = .t.
+  write_Na00  = .t.
+  write_Napj  = .t.
+  write_dens  = .t.
+  write_fvel  = .t.
+  write_temp  = .t.
+/
+&MODEL_PAR
+  ! Collisionality
+  CLOS    = 0
+  NL_CLOS = 0
+  LINEARITY = 'nonlinear'
+  Na      = 1 ! number of species
+  mu_x    = 1.0
+  mu_y    = 1.0
+  N_HD    = 4
+  mu_z    = 2.0
+  HYP_V   = 'hypcoll'
+  mu_p    = 0.0
+  mu_j    = 0.0
+  nu      = 0.001
+  beta    = 0.0
+  ADIAB_E = .t.
+  tau_e   = 1.0
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_  = 1.0
+ sigma_= 1.0
+ q_    = 1.0
+ k_N_  = 2.22
+ k_T_  = 6.96
+/
+
+&COLLISION_PAR
+  collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+  GK_CO           = .f.
+  INTERSPECIES    = .true.
+  mat_file       = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+  INIT_OPT         = 'blob'
+  ACT_ON_MODES     = 'donothing'
+  init_background  = 0.0
+  init_noiselvl    = 0.005
+  iseed            = 42
+/
+&TIME_INTEGRATION_PAR
+  numerical_scheme = 'RK4'
+/
diff --git a/testcases/cyclone_example/fort_01.90 b/testcases/cyclone_example/fort_01.90
new file mode 100644
index 0000000000000000000000000000000000000000..5d11c5ae8792a9a7744030d0ec50bba5522c8110
--- /dev/null
+++ b/testcases/cyclone_example/fort_01.90
@@ -0,0 +1,92 @@
+&BASIC
+  nrun       = 99999999
+  dt         = 0.01
+  tmax       = 75
+  maxruntime = 356400
+  job2load   = 0
+/
+&GRID
+  pmax   = 4
+  jmax   = 2
+  Nx     = 128
+  Lx     = 120
+  Ny     = 48
+  Ly     = 120
+  Nz     = 6!16
+  SG     = .f.
+  Nexc   = 0
+/
+&GEOMETRY
+  geom   = 's-alpha'
+  q0     = 1.4
+  shear  = 0.8
+  eps    = 0.18
+  kappa  = 1.0
+  s_kappa= 0.0
+  delta  = 0.0
+  s_delta= 0.0
+  zeta   = 0.0
+  s_zeta = 0.0
+  parallel_bc = 'dirichlet'
+  shift_y= 0.0
+/
+&OUTPUT_PAR
+  dtsave_0d = 0.1
+  dtsave_1d = -1
+  dtsave_2d = -1
+  dtsave_3d = 1
+  dtsave_5d = 10
+  write_doubleprecision = .f.
+  write_gamma = .t.
+  write_hf    = .t.
+  write_phi   = .t.
+  write_Na00  = .t.
+  write_Napj  = .t.
+  write_dens  = .t.
+  write_fvel  = .t.
+  write_temp  = .t.
+/
+&MODEL_PAR
+  ! Collisionality
+  CLOS    = 0
+  NL_CLOS = 0
+  LINEARITY = 'nonlinear'
+  Na      = 1 ! number of species
+  mu_x    = 1.0
+  mu_y    = 1.0
+  N_HD    = 4
+  mu_z    = 2.0
+  HYP_V   = 'hypcoll'
+  mu_p    = 0.0
+  mu_j    = 0.0
+  nu      = 0.001
+  beta    = 0.0
+  ADIAB_E = .t.
+  tau_e   = 1.0
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_  = 1.0
+ sigma_= 1.0
+ q_    = 1.0
+ k_N_  = 2.22
+ k_T_  = 6.96
+/
+
+&COLLISION_PAR
+  collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+  GK_CO           = .f.
+  INTERSPECIES    = .true.
+  mat_file       = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+  INIT_OPT         = 'blob'
+  ACT_ON_MODES     = 'donothing'
+  init_background  = 0.0
+  init_noiselvl    = 0.005
+  iseed            = 42
+/
+&TIME_INTEGRATION_PAR
+  numerical_scheme = 'RK4'
+/
diff --git a/testcases/cyclone_example/fort_02.90 b/testcases/cyclone_example/fort_02.90
new file mode 100644
index 0000000000000000000000000000000000000000..71f6ce37d2e98c7c59692dce468c95335f19ca1b
--- /dev/null
+++ b/testcases/cyclone_example/fort_02.90
@@ -0,0 +1,92 @@
+&BASIC
+  nrun       = 10000000
+  dt         = 0.01
+  tmax       = 150
+  maxruntime = 356400
+  job2load   = 1
+/
+&GRID
+  pmax   = 4
+  jmax   = 2
+  Nx     = 128
+  Lx     = 120
+  Ny     = 48
+  Ly     = 120
+  Nz     = 16
+  SG     = .f.
+  Nexc   = 0
+/
+&GEOMETRY
+  geom   = 's-alpha'
+  q0     = 1.4
+  shear  = 0.8
+  eps    = 0.18
+  kappa  = 1.0
+  s_kappa= 0.0
+  delta  = 0.0
+  s_delta= 0.0
+  zeta   = 0.0
+  s_zeta = 0.0
+  parallel_bc = 'dirichlet'
+  shift_y= 0.0
+/
+&OUTPUT_PAR
+  dtsave_0d = 1
+  dtsave_1d = -1
+  dtsave_2d = -1
+  dtsave_3d = 1
+  dtsave_5d = 10
+  write_doubleprecision = .f.
+  write_gamma = .t.
+  write_hf    = .t.
+  write_phi   = .t.
+  write_Na00  = .t.
+  write_Napj  = .t.
+  write_dens  = .t.
+  write_fvel  = .t.
+  write_temp  = .t.
+/
+&MODEL_PAR
+  ! Collisionality
+  CLOS    = 0
+  NL_CLOS = 0
+  LINEARITY = 'nonlinear'
+  Na      = 1 ! number of species
+  mu_x    = 1.0
+  mu_y    = 1.0
+  N_HD    = 4
+  mu_z    = 2.0
+  HYP_V   = 'hypcoll'
+  mu_p    = 0.0
+  mu_j    = 0.0
+  nu      = 0.001
+  beta    = 0.0
+  ADIAB_E = .t.
+  tau_e   = 1.0
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_  = 1.0
+ sigma_= 1.0
+ q_    = 1.0
+ k_N_  = 2.22
+ k_T_  = 6.96
+/
+
+&COLLISION_PAR
+  collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+  GK_CO           = .f.
+  INTERSPECIES    = .true.
+  mat_file       = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+  INIT_OPT         = 'blob'
+  ACT_ON_MODES     = 'donothing'
+  init_background  = 0.1
+  init_noiselvl    = 0.00
+  iseed            = 42
+/
+&TIME_INTEGRATION_PAR
+  numerical_scheme = 'RK4'
+/
diff --git a/testcases/cyclone_example/fort_03.90 b/testcases/cyclone_example/fort_03.90
new file mode 100644
index 0000000000000000000000000000000000000000..522faa1ba1fd7bbde75a1e0119dbf30527f744fe
--- /dev/null
+++ b/testcases/cyclone_example/fort_03.90
@@ -0,0 +1,92 @@
+&BASIC
+  nrun       = 10000000
+  dt         = 0.01
+  tmax       = 155
+  maxruntime = 356400
+  job2load   = 2
+/
+&GRID
+  pmax   = 4
+  jmax   = 2
+  Nx     = 128
+  Lx     = 120
+  Ny     = 48
+  Ly     = 120
+  Nz     = 16
+  SG     = .f.
+  Nexc   = 0
+/
+&GEOMETRY
+  geom   = 's-alpha'
+  q0     = 1.4
+  shear  = 0.8
+  eps    = 0.18
+  kappa  = 1.0
+  s_kappa= 0.0
+  delta  = 0.0
+  s_delta= 0.0
+  zeta   = 0.0
+  s_zeta = 0.0
+  parallel_bc = 'dirichlet'
+  shift_y= 0.0
+/
+&OUTPUT_PAR
+  dtsave_0d = 0.01
+  dtsave_1d = -1
+  dtsave_2d = -1
+  dtsave_3d = 1
+  dtsave_5d = 10
+  write_doubleprecision = .f.
+  write_gamma = .t.
+  write_hf    = .t.
+  write_phi   = .t.
+  write_Na00  = .t.
+  write_Napj  = .t.
+  write_dens  = .t.
+  write_fvel  = .t.
+  write_temp  = .t.
+/
+&MODEL_PAR
+  ! Collisionality
+  CLOS    = 0
+  NL_CLOS = 0
+  LINEARITY = 'nonlinear'
+  Na      = 1 ! number of species
+  mu_x    = 1.0
+  mu_y    = 1.0
+  N_HD    = 4
+  mu_z    = 2.0
+  HYP_V   = 'hypcoll'
+  mu_p    = 0.0
+  mu_j    = 0.0
+  nu      = 0.001
+  beta    = 0.0
+  ADIAB_E = .t.
+  tau_e   = 1.0
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_  = 1.0
+ sigma_= 1.0
+ q_    = 1.0
+ k_N_  = 2.22
+ k_T_  = 6.96
+/
+
+&COLLISION_PAR
+  collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+  GK_CO           = .f.
+  INTERSPECIES    = .true.
+  mat_file       = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+  INIT_OPT         = 'blob'
+  ACT_ON_MODES     = 'donothing'
+  init_background  = 0.1
+  init_noiselvl    = 0.00
+  iseed            = 42
+/
+&TIME_INTEGRATION_PAR
+  numerical_scheme = 'RK4'
+/