Added a flag to know if the process has kr=0 data or not

src/grid_mod.F90

@@ -49,8 +49,10 @@ MODULE grid
   REAL(dp), DIMENSION(:),   ALLOCATABLE, PUBLIC :: kparray     ! kperp array
   REAL(dp), PUBLIC, PROTECTED ::  deltakr, deltakz
   INTEGER,  PUBLIC, PROTECTED ::  ikrs, ikre, ikzs, ikze, ikps, ikpe
-  INTEGER,  PUBLIC, PROTECTED :: ikr_0, ikz_0 ! Indices of k-grid origin
+  INTEGER,  PUBLIC, PROTECTED :: ikr_0, ikz_0, ikr_max, ikz_max ! Indices of k-grid origin and max
   INTEGER,  PUBLIC :: ikr, ikz, ip, ij, ikp ! counters
+  LOGICAL,  PUBLIC, PROTECTED :: contains_kr0   = .false. ! rank of the proc containing kr=0 indices
+  LOGICAL,  PUBLIC, PROTECTED :: contains_krmax = .false. ! rank of the proc containing kr=max indices
 
   ! Grid containing the polynomials degrees
   INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_e
@@ -164,10 +166,19 @@ CONTAINS
       deltakr = 2._dp*PI/Lr
     ENDIF
 
-    ! Discretized kr positions ordered as dk*(0 1 2 3)
+    ! Creating a grid ordered as dk*(0 1 2 3)
     DO ikr = ikrs,ikre
       krarray(ikr) = REAL(ikr-1,dp) * deltakr
-      IF (krarray(ikr) .EQ. 0) ikr_0 = ikr
+      ! Finding kr=0
+      IF (krarray(ikr) .EQ. 0) THEN
+        ikr_0 = ikr
+        contains_kr0 = .true.
+      ENDIF
+      ! Finding krmax
+      IF (krarray(ikr) .EQ. (Nr/2+1)*deltakr) THEN
+        ikr_max = ikr
+        contains_krmax = .true.
+      ENDIF
     END DO
 
     ! Orszag 2/3 filter
@@ -193,18 +204,23 @@ CONTAINS
     ikze = Nkz
     ALLOCATE(kzarray(ikzs:ikze))
 
-    ! Grid spacings and discretized kz positions ordered as dk*(0 1 2 3 -2 -1)
     IF (Lz .EQ. 0) THEN ! 1D linear case
       deltakz    = 1._dp
       kzarray(1) = 0
       ikz_0      = 1
+      ikz_max    = 1
     ELSE
       deltakz = 2._dp*PI/Lz
+      ! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
       DO ikz = ikzs,ikze
         kzarray(ikz) = deltakz*(MODULO(ikz-1,Nkz/2)-Nkz/2*FLOOR(2.*real(ikz-1)/real(Nkz)))
         if (ikz .EQ. Nz/2+1)     kzarray(ikz) = -kzarray(ikz)
+        ! Finding kz=0
         IF (kzarray(ikz) .EQ. 0) ikz_0 = ikz
+        ! Finding kzmax
+        IF (kzarray(ikr) .EQ. (Nz/2)*deltakr) ikr_max = ikr
       END DO
+      IF(my_id .EQ. 0) write(*,*) kzarray
     ENDIF
 
     ! Orszag 2/3 filter

src/inital.F90

@@ -100,9 +100,9 @@ SUBROUTINE init_moments
           END DO
         END DO
 
-        IF ( ikrs .EQ. 1 ) THEN
+        IF ( contains_kr0 ) THEN
           DO ikz=2,Nkz/2 !symmetry at kr = 0
-            moments_e( ip,ij,1,ikz, :) = moments_e( ip,ij,1,Nkz+2-ikz, :)
+            moments_e( ip,ij,ikr_0,ikz, :) = moments_e( ip,ij,ikr_0,Nkz+2-ikz, :)
           END DO
         ENDIF
 
@@ -119,9 +119,9 @@ SUBROUTINE init_moments
           END DO
         END DO
 
-        IF ( ikrs .EQ. 1 ) THEN
+        IF ( contains_kr0 ) THEN
           DO ikz=2,Nkz/2 !symmetry at kr = 0
-            moments_i( ip,ij,1,ikz, :) = moments_i( ip,ij,1,Nkz+2-ikz, :)
+            moments_i( ip,ij,ikr_0,ikz, :) = moments_i( ip,ij,ikr_0,Nkz+2-ikz, :)
           END DO
         ENDIF
 
@@ -184,10 +184,11 @@ SUBROUTINE init_phi
       END DO
 
       !symmetry at kr = 0 to keep real inverse transform
-      IF ( ikrs .EQ. 1 ) THEN
+      IF ( contains_kr0 ) THEN
         DO ikz=2,Nkz/2
-          phi(1,ikz) = phi(1,Nkz+2-ikz)
+          phi(ikr_0,ikz) = phi(ikr_0,Nkz+2-ikz)
         END DO
+        phi(ikr_0,Nz/2) = REAL(phi(ikr_0,Nz/2)) !origin must be real
       ENDIF
 
       !**** Cancel previous moments initialization

src/stepon.F90

@@ -129,34 +129,34 @@ SUBROUTINE stepon
       END SUBROUTINE anti_aliasing
 
       SUBROUTINE enforce_symetry ! Force X(k) = X(N-k)* complex conjugate symmetry
-        ! Electron moments
-        DO ip=ips_e,ipe_e
-          DO ij=ijs_e,ije_e
-            DO ikz=2,Nkz/2 !symmetry at kr = 0
-              moments_e( ip,ij,1,ikz, :) = CONJG(moments_e( ip,ij,1,Nkz+2-ikz, :))
+        IF ( contains_kr0 ) THEN
+          ! Electron moments
+          DO ip=ips_e,ipe_e
+            DO ij=ijs_e,ije_e
+              DO ikz=2,Nkz/2 !symmetry at kr = 0
+                moments_e( ip,ij,ikr_0,ikz, :) = CONJG(moments_e( ip,ij,ikr_0,Nkz+2-ikz, :))
+              END DO
+              ! must be real at origin
+              moments_e(ip,ij, ikr_0,ikz_0, :) = REAL(moments_e(ip,ij, ikr_0,ikz_0, :))
             END DO
-            ! must be real at origin and top right
-            moments_e(ip,ij, ikr_0,ikz_0, :) = REAL(moments_e(ip,ij, ikr_0,ikz_0, :))
-            moments_e(ip,ij,  Nr/2,Nz/2 , :) = REAL(moments_e(ip,ij,  Nr/2,Nz/2 , :))
           END DO
-        END DO
-        ! Ion moments
-        DO ip=ips_i,ipe_i
-          DO ij=ijs_i,ije_i
-            DO ikz=2,Nkz/2 !symmetry at kr = 0
-              moments_i( ip,ij,1,ikz, :) = CONJG(moments_i( ip,ij,1,Nkz+2-ikz, :))
+          ! Ion moments
+          DO ip=ips_i,ipe_i
+            DO ij=ijs_i,ije_i
+              DO ikz=2,Nkz/2 !symmetry at kr = 0
+                moments_i( ip,ij,ikr_0,ikz, :) = CONJG(moments_i( ip,ij,ikr_0,Nkz+2-ikz, :))
+              END DO
+              ! must be real at origin and top right
+              moments_i(ip,ij, ikr_0,ikz_0, :) = REAL(moments_i(ip,ij, ikr_0,ikz_0, :))
             END DO
-            ! must be real at origin and top right
-            moments_i(ip,ij, ikr_0,ikz_0, :) = REAL(moments_i(ip,ij, ikr_0,ikz_0, :))
-            moments_i(ip,ij,  Nr/2,Nz/2 , :) = REAL(moments_i(ip,ij,  Nr/2,Nz/2 , :))
           END DO
-        END DO
-        ! Phi
-        DO ikz=2,Nkz/2 !symmetry at kr = 0
-          phi(1,ikz) = phi(1,Nkz+2-ikz)
-        END DO
-        ! must be real at origin and top right
-        phi(ikr_0,ikz_0) = REAL(phi(ikr_0,ikz_0))
+          ! Phi
+          DO ikz=2,Nkz/2 !symmetry at kr = 0
+            phi(ikr_0,ikz) = phi(ikr_0,Nkz+2-ikz)
+          END DO
+          ! must be real at origin
+          phi(ikr_0,ikz_0) = REAL(phi(ikr_0,ikz_0))
+        ENDIF
       END SUBROUTINE enforce_symetry
 
 END SUBROUTINE stepon