diff --git a/src/inital.F90 b/src/inital.F90
index 2dfd084e44a606a3e5d898bc79ccb6229d5c277f..a2273ff64f7c5d9cb5930d53bfd5acf0b2d02341 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -13,6 +13,7 @@ SUBROUTINE inital
USE closure
USE ghosts
USE restarts
+ USE numerics, ONLY: wipe_turbulence, wipe_zonalflow
IMPLICIT NONE
CALL set_updatetlevel(1)
@@ -27,25 +28,33 @@ SUBROUTINE inital
ELSE
! set phi with noise and set moments to 0
IF (INIT_NOISY_PHI) THEN
+ IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
CALL init_phi
-
! set moments_00 (GC density) with noise and compute phi afterwards
ELSE
- IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy moments'
- CALL init_moments ! init noisy N_0
+ IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy gyrocenter density'
+ CALL init_gyrodens ! init only gyrocenter density
+ ! CALL init_moments ! init all moments randomly (unadvised)
CALL poisson ! get phi_0 = phi(N_0)
ENDIF
ENDIF
+ ! Option for wiping the ZF modes (ky==0)
+ IF ( WIPE_ZF .GT. 0 ) THEN
+ IF (my_id .EQ. 0) WRITE(*,*) '-Wiping ZF modes'
+ CALL wipe_zonalflow
+ ENDIF
+
! Option for wiping the turbulence and check growth of secondary inst.
- IF ( WIPE_TURB ) THEN
+ IF ( WIPE_TURB .GT. 0 ) THEN
IF (my_id .EQ. 0) WRITE(*,*) '-Wiping turbulence'
CALL wipe_turbulence
ENDIF
+
! Option for initializing a gaussian blob on the zonal profile
IF ( INIT_BLOB ) THEN
IF (my_id .EQ. 0) WRITE(*,*) '--init a blob'
- CALL put_blob
+ CALL initialize_blob
ENDIF
IF (my_id .EQ. 0) WRITE(*,*) 'Apply closure'
@@ -70,7 +79,7 @@ END SUBROUTINE inital
!******************************************************************************!
!******************************************************************************!
-!!!!!!! Initialize the moments randomly
+!!!!!!! Initialize all the moments randomly
!******************************************************************************!
SUBROUTINE init_moments
USE basic
@@ -155,106 +164,162 @@ SUBROUTINE init_moments
END SUBROUTINE init_moments
!******************************************************************************!
-
!******************************************************************************!
-!!!!!!! Initialize a noisy ES potential and cancel the moments
+!!!!!!! Initialize the gyrocenter density randomly
!******************************************************************************!
-SUBROUTINE init_phi
+SUBROUTINE init_gyrodens
USE basic
USE grid
USE fields
USE prec_const
+ USE utility, ONLY: checkfield
USE initial_par
+ USE model, ONLY : NON_LIN
IMPLICIT NONE
REAL(dp) :: noise
REAL(dp) :: kx, ky, sigma, gain, ky_shift
INTEGER, DIMENSION(12) :: iseedarr
- IF (INIT_NOISY_PHI) THEN
- IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
- ! Seed random number generator
- iseedarr(:)=iseed
- CALL RANDOM_SEED(PUT=iseedarr+my_id)
+ ! Seed random number generator
+ iseedarr(:)=iseed
+ CALL RANDOM_SEED(PUT=iseedarr+my_id)
- !**** noise initialization *******************************************
+ !**** Broad noise initialization *******************************************
+ DO ip=ips_e,ipe_e
+ DO ij=ijs_e,ije_e
- DO ikx=ikxs,ikxe
- DO iky=ikys,ikye
- DO iz=izs,ize
- CALL RANDOM_NUMBER(noise)
- phi(ikx,iky,iz) = (init_background + init_noiselvl*(noise-0.5_dp))*AA_x(ikx)*AA_y(iky)
- ENDDO
+ DO ikx=ikxs,ikxe
+ DO iky=ikys,ikye
+ DO iz=izs,ize
+ CALL RANDOM_NUMBER(noise)
+ IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
+ moments_e(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
+ ELSE
+ moments_e(ip,ij,ikx,iky,iz,:) = 0._dp
+ ENDIF
+ END DO
+ END DO
END DO
+
+ IF ( contains_kx0 ) THEN
+ DO iky=2,Nky/2 !symmetry at kx = 0 for all z
+ moments_e(ip,ij,ikx_0,iky,:,:) = moments_e( ip,ij,ikx_0,Nky+2-iky,:, :)
+ END DO
+ ENDIF
+
END DO
+ END DO
- !symmetry at kx = 0 to keep real inverse transform
- IF ( contains_kx0 ) THEN
- DO iky=2,Nky/2
- phi(ikx_0,iky,:) = phi(ikx_0,Nky+2-iky,:)
- END DO
- phi(ikx_0,Ny/2,:) = REAL(phi(ikx_0,Ny/2,:)) !origin must be real
- ENDIF
+ DO ip=ips_i,ipe_i
+ DO ij=ijs_i,ije_i
- !**** ensure no previous moments initialization
- moments_e = 0._dp; moments_i = 0._dp
+ DO ikx=ikxs,ikxe
+ DO iky=ikys,ikye
+ DO iz=izs,ize
+ CALL RANDOM_NUMBER(noise)
+ IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
+ moments_i(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
+ ELSE
+ moments_i(ip,ij,ikx,iky,iz,:) = 0._dp
+ ENDIF
+ END DO
+ END DO
+ END DO
- !**** Zonal Flow initialization *******************************************
- ! put a mode at ikx = mode number + 1, symmetry is already included since kx>=0
- IF(INIT_ZF .GT. 0) THEN
- IF (my_id .EQ. 0) WRITE(*,*) 'Init ZF phi'
- IF( (INIT_ZF+1 .GT. ikxs) .AND. (INIT_ZF+1 .LT. ikxe) ) THEN
- DO iz = izs,ize
- phi(INIT_ZF+1,iky_0,iz) = ZF_AMP*(2._dp*PI)**2/deltakx/deltaky/2._dp * COS((iz-1)/Nz*2._dp*PI)
- moments_i(1,1,INIT_ZF+1,iky_0,iz,:) = kxarray(INIT_ZF+1)**2*phi(INIT_ZF+1,iky_0,iz)* COS((iz-1)/Nz*2._dp*PI)
- moments_e(1,1,INIT_ZF+1,iky_0,iz,:) = 0._dp
- ENDDO
+ IF ( contains_kx0 ) THEN
+ DO iky=2,Nky/2 !symmetry at kx = 0 for all z
+ moments_i( ip,ij,ikx_0,iky,:,:) = moments_i( ip,ij,ikx_0,Nky+2-iky,:,:)
+ END DO
ENDIF
- ENDIF
- ELSE ! we compute phi from noisy moments and poisson
- CALL poisson
- ENDIF
-END SUBROUTINE init_phi
+ END DO
+ END DO
+
+ ! Putting to zero modes that are not in the 2/3 Orszag rule
+ IF (NON_LIN) THEN
+ DO ikx=ikxs,ikxe
+ DO iky=ikys,ikye
+ DO iz=izs,ize
+ DO ip=ips_e,ipe_e
+ DO ij=ijs_e,ije_e
+ moments_e( ip,ij,ikx,iky,iz, :) = moments_e( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+ DO ip=ips_i,ipe_i
+ DO ij=ijs_i,ije_i
+ moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+END SUBROUTINE init_gyrodens
!******************************************************************************!
!******************************************************************************!
-!!!!!!! Remove all ky!=0 modes to conserve only zonal modes in a restart
+!!!!!!! Initialize a noisy ES potential and cancel the moments
!******************************************************************************!
-SUBROUTINE wipe_turbulence
- USE fields
+SUBROUTINE init_phi
+ USE basic
USE grid
+ USE fields
+ USE prec_const
+ USE initial_par
IMPLICIT NONE
- DO ikx=ikxs,ikxe
- DO iky=ikys,ikye
- DO iz=izs,ize
- DO ip=ips_e,ipe_e
- DO ij=ijs_e,ije_e
- IF( iky .NE. iky_0) THEN
- moments_e( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_e( ip,ij,ikx,iky,iz, :)
- ELSE
- moments_e( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_e( ip,ij,ikx,iky,iz, :)
- ENDIF
- ENDDO
- ENDDO
- DO ip=ips_i,ipe_i
- DO ij=ijs_i,ije_i
- IF( iky .NE. iky_0) THEN
- moments_i( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_i( ip,ij,ikx,iky,iz, :)
- ELSE
- moments_i( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_i( ip,ij,ikx,iky,iz, :)
+
+ REAL(dp) :: noise
+ REAL(dp) :: kx, ky, sigma, gain, ky_shift
+ INTEGER, DIMENSION(12) :: iseedarr
+
+ ! Seed random number generator
+ iseedarr(:)=iseed
+ CALL RANDOM_SEED(PUT=iseedarr+my_id)
+
+ !**** noise initialization *******************************************
+
+ DO ikx=ikxs,ikxe
+ DO iky=ikys,ikye
+ DO iz=izs,ize
+ CALL RANDOM_NUMBER(noise)
+ phi(ikx,iky,iz) = (init_background + init_noiselvl*(noise-0.5_dp))!*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ END DO
+ END DO
+
+ !symmetry at kx = 0 to keep real inverse transform
+ IF ( contains_kx0 ) THEN
+ DO iky=2,Nky/2
+ phi(ikx_0,iky,:) = phi(ikx_0,Nky+2-iky,:)
+ END DO
+ phi(ikx_0,Ny/2,:) = REAL(phi(ikx_0,Ny/2,:)) !origin must be real
+ ENDIF
+
+ !**** ensure no previous moments initialization
+ moments_e = 0._dp; moments_i = 0._dp
+
+ !**** Zonal Flow initialization *******************************************
+ ! put a mode at ikx = mode number + 1, symmetry is already included since kx>=0
+ IF(INIT_ZF .GT. 0) THEN
+ IF (my_id .EQ. 0) WRITE(*,*) 'Init ZF phi'
+ IF( (INIT_ZF+1 .GT. ikxs) .AND. (INIT_ZF+1 .LT. ikxe) ) THEN
+ DO iz = izs,ize
+ phi(INIT_ZF+1,iky_0,iz) = ZF_AMP*(2._dp*PI)**2/deltakx/deltaky/2._dp * COS((iz-1)/Nz*2._dp*PI)
+ moments_i(1,1,INIT_ZF+1,iky_0,iz,:) = kxarray(INIT_ZF+1)**2*phi(INIT_ZF+1,iky_0,iz)* COS((iz-1)/Nz*2._dp*PI)
+ moments_e(1,1,INIT_ZF+1,iky_0,iz,:) = 0._dp
+ ENDDO
ENDIF
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
-END SUBROUTINE
+ ENDIF
+
+END SUBROUTINE init_phi
+!******************************************************************************!
+
!******************************************************************************!
!******************************************************************************!
!!!!!!! Initialize an ionic Gaussian blob on top of the preexisting modes
!******************************************************************************!
-SUBROUTINE put_blob
+SUBROUTINE initialize_blob
USE fields
USE grid
USE model, ONLY: sigmai2_taui_o2
@@ -281,5 +346,5 @@ SUBROUTINE put_blob
ENDDO
ENDDO
ENDDO
-END SUBROUTINE put_blob
+END SUBROUTINE initialize_blob
!******************************************************************************!
diff --git a/src/initial_par_mod.F90 b/src/initial_par_mod.F90
index fefa91183b3087a6c0d8d90bd09f4e4c0238a867..de90b696b7ad4bbd51dfc58b409195049060b78d 100644
--- a/src/initial_par_mod.F90
+++ b/src/initial_par_mod.F90
@@ -10,8 +10,10 @@ MODULE initial_par
! Initialization through a zonal flow phi
INTEGER, PUBLIC, PROTECTED :: INIT_ZF = 0
REAL(DP), PUBLIC, PROTECTED :: ZF_AMP = 1E+3_dp
- ! Wipe turbulence in the restart
- LOGICAL, PUBLIC, PROTECTED :: WIPE_TURB = .false.
+ ! Wipe zonal flow in the restart (=1) or at each step (=2)
+ INTEGER, PUBLIC, PROTECTED :: WIPE_ZF = 0
+ ! Wipe turbulence in the restart (=1) or at each step (=2)
+ INTEGER, PUBLIC, PROTECTED :: WIPE_TURB = 0
! Init a Gaussian blob density in the middle
LOGICAL, PUBLIC, PROTECTED :: INIT_BLOB = .false.
! Initial background level
@@ -37,6 +39,7 @@ CONTAINS
NAMELIST /INITIAL_CON/ INIT_NOISY_PHI
NAMELIST /INITIAL_CON/ INIT_ZF
+ NAMELIST /INITIAL_CON/ WIPE_ZF
NAMELIST /INITIAL_CON/ WIPE_TURB
NAMELIST /INITIAL_CON/ INIT_BLOB
NAMELIST /INITIAL_CON/ init_background
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 1289550f742959fb3ae8d60ecc85a4b298324ef9..e8f1b48692334c667a158e28cb6330d896233944 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -7,6 +7,7 @@ MODULE numerics
implicit none
PUBLIC :: compute_derivatives, build_dnjs_table, evaluate_kernels, compute_lin_coeff
+ PUBLIC :: wipe_turbulence, wipe_zonalflow
CONTAINS
@@ -252,5 +253,79 @@ SUBROUTINE compute_lin_coeff
END SUBROUTINE compute_lin_coeff
+!******************************************************************************!
+!!!!!!! Remove all ky!=0 modes to conserve only zonal modes in a restart
+!******************************************************************************!
+SUBROUTINE wipe_turbulence
+ USE fields
+ USE grid
+ IMPLICIT NONE
+ DO ikx=ikxs,ikxe
+ DO iky=ikys,ikye
+ DO iz=izs,ize
+ DO ip=ips_e,ipe_e
+ DO ij=ijs_e,ije_e
+ IF( iky .NE. iky_0) THEN
+ moments_e( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_e( ip,ij,ikx,iky,iz, :)
+ ELSE
+ moments_e( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_e( ip,ij,ikx,iky,iz, :)
+ ENDIF
+ ENDDO
+ ENDDO
+ DO ip=ips_i,ipe_i
+ DO ij=ijs_i,ije_i
+ IF( iky .NE. iky_0) THEN
+ moments_i( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_i( ip,ij,ikx,iky,iz, :)
+ ELSE
+ moments_i( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_i( ip,ij,ikx,iky,iz, :)
+ ENDIF
+ ENDDO
+ ENDDO
+ IF( iky .NE. iky_0) THEN
+ phi(ikx,iky,iz) = 0e-3_dp*phi(ikx,iky,iz)
+ ELSE
+ phi(ikx,iky,iz) = 1e+0_dp*phi(ikx,iky,iz)
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+END SUBROUTINE
+!******************************************************************************!
+!!!!!!! Remove all ky==0 modes to conserve only non zonal modes in a restart
+!******************************************************************************!
+SUBROUTINE wipe_zonalflow
+ USE fields
+ USE grid
+ IMPLICIT NONE
+ DO ikx=ikxs,ikxe
+ DO iky=ikys,ikye
+ DO iz=izs,ize
+ DO ip=ips_e,ipe_e
+ DO ij=ijs_e,ije_e
+ IF( iky .EQ. iky_0) THEN
+ moments_e( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_e( ip,ij,ikx,iky,iz, :)
+ ELSE
+ moments_e( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_e( ip,ij,ikx,iky,iz, :)
+ ENDIF
+ ENDDO
+ ENDDO
+ DO ip=ips_i,ipe_i
+ DO ij=ijs_i,ije_i
+ IF( iky .EQ. iky_0) THEN
+ moments_i( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_i( ip,ij,ikx,iky,iz, :)
+ ELSE
+ moments_i( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_i( ip,ij,ikx,iky,iz, :)
+ ENDIF
+ ENDDO
+ ENDDO
+ IF( iky .EQ. iky_0) THEN
+ phi(ikx,iky,iz) = 0e-3_dp*phi(ikx,iky,iz)
+ ELSE
+ phi(ikx,iky,iz) = 1e+0_dp*phi(ikx,iky,iz)
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+END SUBROUTINE
END MODULE numerics
diff --git a/src/stepon.F90 b/src/stepon.F90
index 882e1277ccfde6872cdd3f8a413ca48c09b503d5..b7a36785f7d98672baa576a209aad2a674d2f970 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -6,11 +6,13 @@ SUBROUTINE stepon
USE closure
USE collision, ONLY : compute_TColl
USE fields, ONLY: moments_e, moments_i, phi
+ USE initial_par, ONLY: WIPE_ZF, WIPE_TURB
USE ghosts
USE grid
USE model
use prec_const
USE time_integration
+ USE numerics, ONLY: wipe_zonalflow, wipe_turbulence
USE utility, ONLY: checkfield
IMPLICIT NONE
@@ -39,8 +41,11 @@ SUBROUTINE stepon
! Update electrostatic potential phi_n = phi(N_n+1)
CALL poisson
! Update nonlinear term S_n -> S_n+1(phi_n+1,N_n+1)
- IF ( NON_LIN ) CALL compute_Sapj
-
+ IF ( NON_LIN ) CALL compute_Sapj
+ ! Cancel zonal modes artificially
+ IF ( WIPE_ZF .EQ. 2) CALL wipe_zonalflow
+ ! Cancel non zonal modes artificially
+ IF ( WIPE_TURB .EQ. 2) CALL wipe_turbulence
!- Check before next step
CALL checkfield_all()
IF( nlend ) EXIT ! exit do loop