From 7f510053175d3f2d6342eb63d77087b13d7abd09 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Mon, 7 Nov 2022 15:45:42 +0100
Subject: [PATCH] removing unused variables and archiving routines in a new
 source file

---
 src/numerical_experiments_mod.F90 | 113 ++++++++++++++++++++++++++++
 src/numerics_mod.F90              | 120 +++++-------------------------
 2 files changed, 131 insertions(+), 102 deletions(-)
 create mode 100644 src/numerical_experiments_mod.F90

diff --git a/src/numerical_experiments_mod.F90 b/src/numerical_experiments_mod.F90
new file mode 100644
index 00000000..8cdefd0f
--- /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 4c11273e..09fed922 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
-- 
GitLab