From ba6c6cad36e972aae32bb57e510571a4ed3a7ce5 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 3 Feb 2023 16:26:32 +0100
Subject: [PATCH] esthetic and Napjz are real arrays now

---
 src/array_mod.F90          |  4 ++--
 src/diagnose.F90           | 20 ++++++-----------
 src/moments_eq_rhs_mod.F90 | 32 +++++++++++++-------------
 src/parallel_mod.F90       | 46 ++++++++++++++++++--------------------
 src/processing_mod.F90     | 12 +++++-----
 5 files changed, 54 insertions(+), 60 deletions(-)

diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 56fc14b9..4f4096df 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -73,8 +73,8 @@ MODULE array
   COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ni00
 
   ! Kinetic spectrum sum_kx,ky(|Napj(z)|^2), (ip,ij,iz) (should be real)
-  COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Nepjz
-  COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Nipjz
+  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: Nepjz
+  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: Nipjz
 
   ! particle density for electron and ions (iky,ikx,iz)
   COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: dens_e
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index b3e53fc2..1863dccf 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -1,6 +1,5 @@
 SUBROUTINE diagnose(kstep)
   !   Diagnostics, writing simulation state to disk
-
   USE basic
   USE diagnostics_par
   USE processing, ONLY: gflux_ri, hflux_xi
@@ -29,8 +28,6 @@ SUBROUTINE diagnose(kstep)
   CALL diagnose_full(kstep)
 
   CALL cpu_time(t1_diag); tc_diag = tc_diag + (t1_diag - t0_diag)
-
-
 END SUBROUTINE diagnose
 
 SUBROUTINE init_outfile(comm,file0,file,fid)
@@ -333,7 +330,6 @@ SUBROUTINE diagnose_full(kstep)
      CALL closef(fidres)
 
   END IF
-
 END SUBROUTINE diagnose_full
 
 !!-------------- Auxiliary routines -----------------!!
@@ -388,7 +384,6 @@ SUBROUTINE diagnose_0d
 END SUBROUTINE diagnose_0d
 
 SUBROUTINE diagnose_3d
-
   USE basic
   USE futils, ONLY: append, getatt, attach, putarrnd, putarr
   USE fields
@@ -399,7 +394,6 @@ SUBROUTINE diagnose_3d
   USE prec_const
   USE processing
   USE model, ONLY: KIN_E
-
   IMPLICIT NONE
 
   CALL append(fidres,  "/data/var3d/time",           time,ionode=0)
@@ -484,16 +478,16 @@ SUBROUTINE diagnose_3d
   SUBROUTINE write_field3d_pjz_i(field, text)
     USE parallel, ONLY : gather_pjz_i
     IMPLICIT NONE
-    COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize), INTENT(IN) :: field
-    COMPLEX(dp), DIMENSION(1:pmaxi+1,1:jmaxi+1,1:Nz) :: field_full
+    REAL(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize), INTENT(IN) :: field
+    REAL(dp), DIMENSION(1:Np_i,1:Nj_i,1:Nz) :: field_full
     CHARACTER(*), INTENT(IN) :: text
     CHARACTER(LEN=50) :: dset_name
     WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
     IF (num_procs .EQ. 1) THEN ! no data distribution
       CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,izs:ize), ionode=0)
     ELSE
-      CALL gather_pjz_i(field(ips_i:ipe_i,ijs_i:ije_i,izs:ize),field_full(1:pmaxi+1,1:jmaxi+1,1:Nz))
-      CALL putarr(fidres, dset_name, field(1:pmaxi+1,1:jmaxi+1,1:Nz), ionode=0)
+      CALL gather_pjz_i(field(ips_i:ipe_i,ijs_i:ije_i,izs:ize),field_full(1:Np_i,1:Nj_i,1:Nz))
+      CALL putarr(fidres, dset_name, field(1:Np_i,1:Nj_i,1:Nz), ionode=0)
     ENDIF
     CALL attach(fidres, dset_name, "time", time)
   END SUBROUTINE write_field3d_pjz_i
@@ -501,8 +495,8 @@ SUBROUTINE diagnose_3d
   SUBROUTINE write_field3d_pjz_e(field, text)
     USE parallel, ONLY : gather_pjz_e
     IMPLICIT NONE
-    COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize), INTENT(IN) :: field
-    COMPLEX(dp), DIMENSION(1:pmaxe+1,1:jmaxe+1,1:Nz) :: field_full
+    REAL(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize), INTENT(IN) :: field
+    REAL(dp), DIMENSION(1:pmaxe+1,1:jmaxe+1,1:Nz) :: field_full
     CHARACTER(*), INTENT(IN) :: text
     CHARACTER(LEN=50) :: dset_name
     WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
@@ -510,7 +504,7 @@ SUBROUTINE diagnose_3d
       CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,izs:ize), ionode=0)
     ELSE
       CALL gather_pjz_e(field(ips_e:ipe_e,ijs_e:ije_e,izs:ize),field_full(1:pmaxe+1,1:jmaxe+1,1:Nz))
-      CALL putarr(fidres, dset_name, field(1:pmaxi+1,1:jmaxi+1,1:Nz), ionode=0)
+      CALL putarr(fidres, dset_name, field(1:Np_i,1:Nj_i,1:Nz), ionode=0)
     ENDIF
     CALL attach(fidres, dset_name, "time", time)
   END SUBROUTINE write_field3d_pjz_e
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index aad64e46..84dbced3 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -66,24 +66,26 @@ SUBROUTINE compute_moments_eq_rhs
     INTEGER, INTENT(IN) :: ips, ipe, ipgs, ipge, ijs, ije, ijgs, ijge
     INTEGER,  DIMENSION(ips:ipe), INTENT(IN) :: parray
     INTEGER,  DIMENSION(ijs:ije), INTENT(IN) :: jarray
-    REAL(dp), DIMENSION(ips:ipe,ijs:ije), INTENT(IN) :: xnapj,&
-                                        ynapp1j, ynapm1j,   ynapp1jm1, ynapm1jm1,&
-                                        znapm1j, znapm1jp1, znapm1jm1,&
-                                        xphij, xphijp1, xphijm1,&
-                                        xpsij, xpsijp1, xpsijm1
-    REAL(dp), DIMENSION(ips:ipe), INTENT(IN) :: &
-                                        xnapp1j, xnapm1j,   xnapp2j,   xnapm2j
-
-    REAL(dp), DIMENSION(ijs:ije), INTENT(IN) :: xnapjp1, xnapjm1
+    REAL(dp), DIMENSION(ips:ipe,ijs:ije), INTENT(IN) :: xnapj
+    REAL(dp), DIMENSION(ips:ipe),         INTENT(IN) :: xnapp2j, xnapm2j
+    REAL(dp), DIMENSION(ijs:ije),         INTENT(IN) :: xnapjp1, xnapjm1
+    REAL(dp), DIMENSION(ips:ipe),         INTENT(IN) :: xnapp1j, xnapm1j
+    REAL(dp), DIMENSION(ips:ipe,ijs:ije), INTENT(IN) :: ynapp1j, ynapp1jm1, ynapm1j, ynapm1jm1
+    REAL(dp), DIMENSION(ips:ipe,ijs:ije), INTENT(IN) :: znapm1j, znapm1jp1, znapm1jm1
+    REAL(dp), DIMENSION(ips:ipe,ijs:ije), INTENT(IN) :: xphij, xphijp1, xphijm1
+    REAL(dp), DIMENSION(ips:ipe,ijs:ije), INTENT(IN) :: xpsij, xpsijp1, xpsijm1
 
     REAL(dp), DIMENSION(ijgs:ijge,ikys:ikye,ikxs:ikxe,izgs:izge,0:1),INTENT(IN) :: kernel
 
-    COMPLEX(dp), DIMENSION(ipgs:ipge,ijgs:ijge,ikys:ikye,ikxs:ikxe,izgs:izge),INTENT(IN) ::&
-        moments_,nadiab_moments, ddz_napj, interp_napj, ddzND_napj
-    COMPLEX(dp), DIMENSION(ips:ipe,ijs:ije,ikys:ikye,ikxs:ikxe,izs:ize),INTENT(IN) :: Sapj, TColl_
-
+    COMPLEX(dp), DIMENSION(ipgs:ipge,ijgs:ijge,ikys:ikye,ikxs:ikxe,izgs:izge),INTENT(IN) :: nadiab_moments
+    COMPLEX(dp), DIMENSION(ipgs:ipge,ijgs:ijge,ikys:ikye,ikxs:ikxe,izgs:izge),INTENT(IN) :: ddz_napj
+    COMPLEX(dp), DIMENSION(ipgs:ipge,ijgs:ijge,ikys:ikye,ikxs:ikxe,izgs:izge),INTENT(IN) :: interp_napj
+    COMPLEX(dp), DIMENSION( ips:ipe,  ijs:ije, ikys:ikye,ikxs:ikxe, izs:ize), INTENT(IN) :: Sapj
+    COMPLEX(dp), DIMENSION(ipgs:ipge,ijgs:ijge,ikys:ikye,ikxs:ikxe,izgs:izge),INTENT(IN) :: moments_
+    COMPLEX(dp), DIMENSION( ips:ipe,  ijs:ije, ikys:ikye,ikxs:ikxe, izs:ize), INTENT(IN) :: TColl_
+    COMPLEX(dp), DIMENSION(ipgs:ipge,ijgs:ijge,ikys:ikye,ikxs:ikxe,izgs:izge),INTENT(IN) :: ddzND_napj
     !! OUTPUT
-    COMPLEX(dp), DIMENSION(ips:ipe,ijs:ije,ikys:ikye,ikxs:ikxe,izs:ize),INTENT(OUT) :: moments_rhs_
+    COMPLEX(dp), DIMENSION( ips:ipe,  ijs:ije, ikys:ikye,ikxs:ikxe, izs:ize),INTENT(OUT) :: moments_rhs_
 
     INTEGER     :: p_int, j_int ! loops indices and polynom. degrees
     REAL(dp)    :: kx, ky, kperp2
@@ -117,7 +119,7 @@ SUBROUTINE compute_moments_eq_rhs
               eo    = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
               kperp2= kparray(iky,ikx,iz,eo)**2
 
-            IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN
+            IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN ! for the closure scheme
               !! Compute moments_ mixing terms
               Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp
               ! Perpendicular dynamic
diff --git a/src/parallel_mod.F90 b/src/parallel_mod.F90
index 94138539..aa3ed37a 100644
--- a/src/parallel_mod.F90
+++ b/src/parallel_mod.F90
@@ -25,7 +25,6 @@ CONTAINS
 
   SUBROUTINE init_parallel_var
     INTEGER :: i_
-
     !!!!!! XYZ gather variables
     !! Y reduction at constant x and z
     ! number of points recieved and displacement for the y reduction
@@ -122,7 +121,6 @@ CONTAINS
     IF(KIN_E) &
     ALLOCATE(buff_pjxy_zBC_e(ipgs_e:ipge_e,ijgs_e:ijge_e,ikys:ikye,ikxs:ikxe,-2:2))
     ALLOCATE(buff_pjxy_zBC_i(ipgs_i:ipge_i,ijgs_i:ijge_i,ikys:ikye,ikxs:ikxe,-2:2))
-
   END SUBROUTINE init_parallel_var
 
   !!!!! Gather a field in spatial coordinates on rank 0 !!!!!
@@ -165,12 +163,12 @@ CONTAINS
 
   !!!!! Gather a field in kinetic + z coordinates on rank 0 !!!!!
   SUBROUTINE gather_pjz_i(field_sub,field_full)
-    COMPLEX(dp), DIMENSION(ips_i:ipe_i,   1:jmaxi+1, izs:ize), INTENT(IN)    :: field_sub
-    COMPLEX(dp), DIMENSION(   1:pmaxi+1,  1:jmaxi+1,   1:Nz),  INTENT(INOUT) :: field_full
-    COMPLEX(dp), DIMENSION(ips_i:ipe_i)             :: buffer_lp_cz !local p, constant z
-    COMPLEX(dp), DIMENSION(   1:pmaxi+1 )           :: buffer_fp_cz !full  p, constant z
-    COMPLEX(dp), DIMENSION(   1:pmaxi+1,  izs:ize ) :: buffer_fp_lz !full  p, local z
-    COMPLEX(dp), DIMENSION(   1:pmaxi+1,    1:Nz  ) :: buffer_fp_fz !full  p, full  z
+    REAL(dp), DIMENSION(ips_i:ipe_i, 1:Nj_i, izs:ize), INTENT(IN)    :: field_sub
+    REAL(dp), DIMENSION(   1:Np_i,   1:Nj_i,   1:Nz),  INTENT(INOUT) :: field_full
+    REAL(dp), DIMENSION(ips_i:ipe_i)          :: buffer_lp_cz !local p, constant z
+    REAL(dp), DIMENSION(   1:Np_i )           :: buffer_fp_cz !full  p, constant z
+    REAL(dp), DIMENSION(   1:Np_i,  izs:ize ) :: buffer_fp_lz !full  p, local z
+    REAL(dp), DIMENSION(   1:Np_i,    1:Nz  ) :: buffer_fp_fz !full  p, full  z
     INTEGER :: snd_p, snd_z, root_p, root_z, root_ky, ij, iz
 
     snd_p  = local_np_i    ! Number of points to send along p (per z)
@@ -178,20 +176,20 @@ CONTAINS
 
     root_p = 0; root_z = 0; root_ky = 0
     IF(rank_ky .EQ. root_ky) THEN
-      DO ij = 1,jmaxi+1
+      DO ij = 1,Nj_i
         DO iz = izs,ize
           ! fill a buffer to contain a slice of data at constant j and z
           buffer_lp_cz(ips_i:ipe_i) = field_sub(ips_i:ipe_i,ij,iz)
-          CALL MPI_GATHERV(buffer_lp_cz, snd_p,            MPI_DOUBLE_COMPLEX, &
-                           buffer_fp_cz, rcv_p_i, dsp_p_i, MPI_DOUBLE_COMPLEX, &
+          CALL MPI_GATHERV(buffer_lp_cz, snd_p,            MPI_DOUBLE_PRECISION, &
+                           buffer_fp_cz, rcv_p_i, dsp_p_i, MPI_DOUBLE_PRECISION, &
                            root_p, comm_p, ierr)
           buffer_fp_lz(1:Np_i,iz) = buffer_fp_cz(1:Np_i)
         ENDDO
 
         ! send the full line on y contained by root_kyas
         IF(rank_p .EQ. 0) THEN
-          CALL MPI_GATHERV(buffer_fp_lz, snd_z,              MPI_DOUBLE_COMPLEX, &
-                           buffer_fp_fz, rcv_zp_i, dsp_zp_i, MPI_DOUBLE_COMPLEX, &
+          CALL MPI_GATHERV(buffer_fp_lz, snd_z,              MPI_DOUBLE_PRECISION, &
+                           buffer_fp_fz, rcv_zp_i, dsp_zp_i, MPI_DOUBLE_PRECISION, &
                            root_z, comm_z, ierr)
         ENDIF
         ! ID 0 (the one who output) rebuild the whole array
@@ -202,12 +200,12 @@ CONTAINS
   END SUBROUTINE gather_pjz_i
 
   SUBROUTINE gather_pjz_e(field_sub,field_full)
-    COMPLEX(dp), DIMENSION(ips_e:ipe_e,   1:jmaxe+1, izs:ize), INTENT(IN)    :: field_sub
-    COMPLEX(dp), DIMENSION(   1:pmaxe+1,  1:jmaxe+1,   1:Nz),  INTENT(INOUT) :: field_full
-    COMPLEX(dp), DIMENSION(ips_e:ipe_e)         :: buffer_lp_cz !local p, constant z
-    COMPLEX(dp), DIMENSION(   1:pmaxe+1 )       :: buffer_fp_cz !full  p, constant z
-    COMPLEX(dp), DIMENSION(   1:pmaxe+1,  izs:ize ) :: buffer_fp_lz !full  p, local z
-    COMPLEX(dp), DIMENSION(   1:pmaxe+1,    1:Nz  ) :: buffer_fp_fz !full  p, full  z
+    REAL(dp), DIMENSION(ips_e:ipe_e,   1:jmaxe+1, izs:ize), INTENT(IN)    :: field_sub
+    REAL(dp), DIMENSION(   1:pmaxe+1,  1:jmaxe+1,   1:Nz),  INTENT(INOUT) :: field_full
+    REAL(dp), DIMENSION(ips_e:ipe_e)         :: buffer_lp_cz !local p, constant z
+    REAL(dp), DIMENSION(   1:pmaxe+1 )       :: buffer_fp_cz !full  p, constant z
+    REAL(dp), DIMENSION(   1:pmaxe+1,  izs:ize ) :: buffer_fp_lz !full  p, local z
+    REAL(dp), DIMENSION(   1:pmaxe+1,    1:Nz  ) :: buffer_fp_fz !full  p, full  z
     INTEGER :: snd_p, snd_z, root_p, root_z, root_ky, ij, iz
 
     snd_p  = local_np_e    ! Number of points to send along p (per z)
@@ -215,20 +213,20 @@ CONTAINS
 
     root_p = 0; root_z = 0; root_ky = 0
     IF(rank_ky .EQ. root_ky) THEN
-      DO ij = 1,jmaxi+1
+      DO ij = 1,Nj_i
         DO iz = izs,ize
           ! fill a buffer to contain a slice of data at constant j and z
           buffer_lp_cz(ips_e:ipe_e) = field_sub(ips_e:ipe_e,ij,iz)
-          CALL MPI_GATHERV(buffer_lp_cz, snd_p,            MPI_DOUBLE_COMPLEX, &
-                           buffer_fp_cz, rcv_p_e, dsp_p_e, MPI_DOUBLE_COMPLEX, &
+          CALL MPI_GATHERV(buffer_lp_cz, snd_p,            MPI_DOUBLE_PRECISION, &
+                           buffer_fp_cz, rcv_p_e, dsp_p_e, MPI_DOUBLE_PRECISION, &
                            root_p, comm_p, ierr)
           buffer_fp_lz(1:Np_e,iz) = buffer_fp_cz(1:Np_e)
         ENDDO
 
         ! send the full line on y contained by root_kyas
         IF(rank_p .EQ. 0) THEN
-          CALL MPI_GATHERV(buffer_fp_lz, snd_z,              MPI_DOUBLE_COMPLEX, &
-                           buffer_fp_fz, rcv_zp_e, dsp_zp_e, MPI_DOUBLE_COMPLEX, &
+          CALL MPI_GATHERV(buffer_fp_lz, snd_z,              MPI_DOUBLE_PRECISION, &
+                           buffer_fp_fz, rcv_zp_e, dsp_zp_e, MPI_DOUBLE_PRECISION, &
                            root_z, comm_z, ierr)
         ENDIF
         ! ID 0 (the one who output) rebuild the whole array
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index aceca824..e8589da9 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -535,8 +535,8 @@ SUBROUTINE compute_Napjz_spectrum
   USE array,  ONLY : Nipjz, Nepjz
   USE time_integration, ONLY : updatetlevel
   IMPLICIT NONE
-  COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize) :: local_sum_e,global_sum_e, buffer_e
-  COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize) :: local_sum_i,global_sum_i, buffer_i
+  REAL(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize) :: local_sum_e,global_sum_e, buffer_e
+  REAL(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize) :: local_sum_i,global_sum_i, buffer_i
   INTEGER  :: i_, root, count
   root = 0
   ! Electron moments spectrum
@@ -546,8 +546,8 @@ SUBROUTINE compute_Napjz_spectrum
     DO ikx = ikxs,ikxe
       DO iky = ikys,ikye
         local_sum_e(ips_e:ipe_e,ijs_e:ije_e,izs:ize)  = local_sum_e(ips_e:ipe_e,ijs_e:ije_e,izs:ize)  + &
-        (moments_e(ips_e:ips_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel)&
-         * CONJG(moments_e(ips_e:ips_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel)))
+        REAL(REAL(moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel)&
+         * CONJG(moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel))))
       ENDDO
     ENDDO
     ! sum reduction
@@ -588,12 +588,12 @@ SUBROUTINE compute_Napjz_spectrum
   IF (num_procs_ky .GT. 1) THEN
       !! Everyone sends its local_sum to root = 0
       IF (rank_ky .NE. root) THEN
-          CALL MPI_SEND(buffer_i, count , MPI_DOUBLE_COMPLEX, root, 1234, comm_ky, ierr)
+          CALL MPI_SEND(buffer_i, count , MPI_DOUBLE_PRECISION, root, 5678, comm_ky, ierr)
       ELSE
           ! Recieve from all the other processes
           DO i_ = 0,num_procs_ky-1
               IF (i_ .NE. rank_ky) &
-                  CALL MPI_RECV(buffer_i, count , MPI_DOUBLE_COMPLEX, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
+                  CALL MPI_RECV(buffer_i, count , MPI_DOUBLE_PRECISION, i_, 5678, comm_ky, MPI_STATUS_IGNORE, ierr)
                   global_sum_i = global_sum_i + buffer_i
           ENDDO
       ENDIF
-- 
GitLab