diff --git a/src/numerical_experiments_mod.F90 b/src/numerical_experiments_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..8cdefd0f8343db8395b05a5756fff4abfcc59808
--- /dev/null
+++ b/src/numerical_experiments_mod.F90
@@ -0,0 +1,113 @@
+!! The numerical_experiments module contains routines to "play" with the fourier
+! modes in order to understand mechanisms. These routines are not integrated in
+! the main code anymore because they are not used. This file serves as an archive.
+MODULE numerical_experiments
+USE basic
+USE prec_const
+USE grid
+USE utility
+
+implicit none
+
+PUBLIC :: play_with_modes, save_EM_ZF_modes
+
+CONTAINS
+!******************************************************************************!
+!!!!!!! Routine that can artificially increase or wipe modes
+!******************************************************************************!
+SUBROUTINE save_EM_ZF_modes
+ USE fields
+ USE array, ONLY : moments_e_ZF, moments_i_ZF, phi_ZF, moments_e_EM,moments_i_EM,phi_EM
+ USE grid
+ USE time_integration, ONLY: updatetlevel
+ USE model, ONLY: KIN_E
+ IMPLICIT NONE
+ ! Store Zonal and entropy modes
+ IF(contains_ky0) THEN
+ IF(KIN_E) &
+ moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel)
+ moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel)
+ phi_ZF(ikxs:ikxe,izs:ize) = phi(iky_0,ikxs:ikxe,izs:ize)
+ ELSE
+ IF(KIN_E) &
+ moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) = 0._dp
+ moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) = 0._dp
+ phi_ZF(ikxs:ikxe,izs:ize) = 0._dp
+ ENDIF
+ IF(contains_kx0) THEN
+ IF(KIN_E) &
+ moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikx_0,izs:ize,updatetlevel)
+ moments_i_EM(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikx_0,izs:ize,updatetlevel)
+ phi_EM(ikys:ikye,izs:ize) = phi(ikys:ikye,ikx_0,izs:ize)
+ ELSE
+ IF(KIN_E) &
+ moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = 0._dp
+ moments_i_EM(ips_e:ipe_e,ijs_i:ije_i,ikys:ikye,izs:ize) = 0._dp
+ phi_EM(ikys:ikye,izs:ize) = 0._dp
+ ENDIF
+END SUBROUTINE
+
+SUBROUTINE play_with_modes
+ USE fields
+ USE array, ONLY : moments_e_ZF, moments_i_ZF, phi_ZF, moments_e_EM,moments_i_EM,phi_EM
+ USE grid
+ USE time_integration, ONLY: updatetlevel
+ USE initial_par, ONLY: ACT_ON_MODES
+ USE model, ONLY: KIN_E
+ IMPLICIT NONE
+ REAL(dp) :: AMP = 1.5_dp
+
+ SELECT CASE(ACT_ON_MODES)
+ CASE('wipe_zonal') ! Errase the zonal flow
+ IF(KIN_E) &
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = 0._dp
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = 0._dp
+ phi(iky_0,ikxs:ikxe,izs:ize) = 0._dp
+ CASE('wipe_entropymode')
+ IF(KIN_E) &
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikx_0,izs:ize,updatetlevel) = 0._dp
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikx_0,izs:ize,updatetlevel) = 0._dp
+ phi(ikys:ikye,ikx_0,izs:ize) = 0._dp
+ CASE('wipe_turbulence')
+ DO ikx = ikxs,ikxe
+ DO iky = ikys, ikye
+ IF ( (ikx .NE. ikx_0) .AND. (iky .NE. iky_0) ) THEN
+ IF(KIN_E) &
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel) = 0._dp
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,iky,ikx,izs:ize,updatetlevel) = 0._dp
+ phi(iky,ikx,izs:ize) = 0._dp
+ ENDIF
+ ENDDO
+ ENDDO
+ CASE('wipe_nonzonal')
+ DO ikx = ikxs,ikxe
+ DO iky = ikys, ikye
+ IF ( (ikx .NE. ikx_0) ) THEN
+ IF(KIN_E) &
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel) = 0._dp
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,iky,ikx,izs:ize,updatetlevel) = 0._dp
+ phi(iky,ikx,izs:ize) = 0._dp
+ ENDIF
+ ENDDO
+ ENDDO
+ CASE('freeze_zonal')
+ IF(KIN_E) &
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize)
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize)
+ phi(iky_0,ikxs:ikxe,izs:ize) = phi_ZF(ikxs:ikxe,izs:ize)
+ CASE('freeze_entropymode')
+ IF(contains_kx0) THEN
+ IF(KIN_E) &
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikx_0,izs:ize,updatetlevel) = moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize)
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikx_0,izs:ize,updatetlevel) = moments_i_EM(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,izs:ize)
+ phi(ikys:ikye,ikx_0,izs:ize) = phi_EM(ikys:ikye,izs:ize)
+ ENDIF
+ CASE('amplify_zonal')
+ IF(KIN_E) &
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = AMP*moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize)
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = AMP*moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize)
+ phi(iky_0,ikxs:ikxe,izs:ize) = AMP*phi_ZF(ikxs:ikxe,izs:ize)
+ END SELECT
+END SUBROUTINE
+
+END MODULE numerical experiments
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 4c11273eda1d0dcdfc592a74f5b37aa1515f562b..09fed922b4ee78b94261b46e409df2f8dcb22f78 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -1,3 +1,6 @@
+!! MODULE NUMERICS
+! The module numerics contains a set of routines that are called only once at
+! the begining of a run. These routines do not need to be optimzed
MODULE numerics
USE basic
USE prec_const
@@ -7,7 +10,7 @@ MODULE numerics
implicit none
PUBLIC :: build_dnjs_table, evaluate_kernels, evaluate_EM_op
- PUBLIC :: compute_lin_coeff, play_with_modes, save_EM_ZF_modes
+ PUBLIC :: compute_lin_coeff
CONTAINS
@@ -54,7 +57,7 @@ SUBROUTINE evaluate_kernels
USE model, ONLY : sigmae2_taue_o2, sigmai2_taui_o2, KIN_E
IMPLICIT NONE
INTEGER :: j_int
- REAL(dp) :: j_dp, y_
+ REAL(dp) :: j_dp, y_, factj
DO eo = 0,1
DO ikx = ikxs,ikxe
@@ -66,7 +69,13 @@ DO iz = izgs,izge
j_int = jarray_e(ij)
j_dp = REAL(j_int,dp)
y_ = sigmae2_taue_o2 * kparray(iky,ikx,iz,eo)**2
- kernel_e(ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
+ IF(j_int .LE. 0) THEN
+ factj = 1._dp
+ kernel_e(ij,iky,ikx,iz,eo) = 0._dp
+ ELSE
+ factj = GAMMA(j_dp+1._dp)
+ kernel_e(ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/factj
+ ENDIF
ENDDO
IF (ijs_e .EQ. 1) &
kernel_e(ijgs_e,iky,ikx,iz,eo) = 0._dp
@@ -76,7 +85,12 @@ DO iz = izgs,izge
j_int = jarray_i(ij)
j_dp = REAL(j_int,dp)
y_ = sigmai2_taui_o2 * kparray(iky,ikx,iz,eo)**2
- kernel_i(ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
+ IF(j_int .LT. 0) THEN
+ kernel_i(ij,iky,ikx,iz,eo) = 0._dp
+ ELSE
+ factj = GAMMA(j_dp+1._dp)
+ kernel_i(ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/factj
+ ENDIF
ENDDO
IF (ijs_i .EQ. 1) &
kernel_i(ijgs_i,iky,ikx,iz,eo) = 0._dp
@@ -378,102 +392,4 @@ SUBROUTINE compute_lin_coeff
ENDDO
END SUBROUTINE compute_lin_coeff
-!******************************************************************************!
-!!!!!!! Routine that can artificially increase or wipe modes
-!******************************************************************************!
-SUBROUTINE save_EM_ZF_modes
- USE fields
- USE array, ONLY : moments_e_ZF, moments_i_ZF, phi_ZF, moments_e_EM,moments_i_EM,phi_EM
- USE grid
- USE time_integration, ONLY: updatetlevel
- USE model, ONLY: KIN_E
- IMPLICIT NONE
- ! Store Zonal and entropy modes
- IF(contains_ky0) THEN
- IF(KIN_E) &
- moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel)
- moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel)
- phi_ZF(ikxs:ikxe,izs:ize) = phi(iky_0,ikxs:ikxe,izs:ize)
- ELSE
- IF(KIN_E) &
- moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) = 0._dp
- moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) = 0._dp
- phi_ZF(ikxs:ikxe,izs:ize) = 0._dp
- ENDIF
- IF(contains_kx0) THEN
- IF(KIN_E) &
- moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikx_0,izs:ize,updatetlevel)
- moments_i_EM(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikx_0,izs:ize,updatetlevel)
- phi_EM(ikys:ikye,izs:ize) = phi(ikys:ikye,ikx_0,izs:ize)
- ELSE
- IF(KIN_E) &
- moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = 0._dp
- moments_i_EM(ips_e:ipe_e,ijs_i:ije_i,ikys:ikye,izs:ize) = 0._dp
- phi_EM(ikys:ikye,izs:ize) = 0._dp
- ENDIF
-END SUBROUTINE
-
-SUBROUTINE play_with_modes
- USE fields
- USE array, ONLY : moments_e_ZF, moments_i_ZF, phi_ZF, moments_e_EM,moments_i_EM,phi_EM
- USE grid
- USE time_integration, ONLY: updatetlevel
- USE initial_par, ONLY: ACT_ON_MODES
- USE model, ONLY: KIN_E
- IMPLICIT NONE
- REAL(dp) :: AMP = 1.5_dp
-
- SELECT CASE(ACT_ON_MODES)
- CASE('wipe_zonal') ! Errase the zonal flow
- IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = 0._dp
- moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = 0._dp
- phi(iky_0,ikxs:ikxe,izs:ize) = 0._dp
- CASE('wipe_entropymode')
- IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikx_0,izs:ize,updatetlevel) = 0._dp
- moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikx_0,izs:ize,updatetlevel) = 0._dp
- phi(ikys:ikye,ikx_0,izs:ize) = 0._dp
- CASE('wipe_turbulence')
- DO ikx = ikxs,ikxe
- DO iky = ikys, ikye
- IF ( (ikx .NE. ikx_0) .AND. (iky .NE. iky_0) ) THEN
- IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel) = 0._dp
- moments_i(ips_i:ipe_i,ijs_i:ije_i,iky,ikx,izs:ize,updatetlevel) = 0._dp
- phi(iky,ikx,izs:ize) = 0._dp
- ENDIF
- ENDDO
- ENDDO
- CASE('wipe_nonzonal')
- DO ikx = ikxs,ikxe
- DO iky = ikys, ikye
- IF ( (ikx .NE. ikx_0) ) THEN
- IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel) = 0._dp
- moments_i(ips_i:ipe_i,ijs_i:ije_i,iky,ikx,izs:ize,updatetlevel) = 0._dp
- phi(iky,ikx,izs:ize) = 0._dp
- ENDIF
- ENDDO
- ENDDO
- CASE('freeze_zonal')
- IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize)
- moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize)
- phi(iky_0,ikxs:ikxe,izs:ize) = phi_ZF(ikxs:ikxe,izs:ize)
- CASE('freeze_entropymode')
- IF(contains_kx0) THEN
- IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikx_0,izs:ize,updatetlevel) = moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize)
- moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikx_0,izs:ize,updatetlevel) = moments_i_EM(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,izs:ize)
- phi(ikys:ikye,ikx_0,izs:ize) = phi_EM(ikys:ikye,izs:ize)
- ENDIF
- CASE('amplify_zonal')
- IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = AMP*moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize)
- moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = AMP*moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize)
- phi(iky_0,ikxs:ikxe,izs:ize) = AMP*phi_ZF(ikxs:ikxe,izs:ize)
- END SELECT
-END SUBROUTINE
-
END MODULE numerics