diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 1fed4a4adba0903b215f27965941b8cf44a9aa31..383f85a537a733c189bfd46440d45f8969586440 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -5,7 +5,8 @@ MODULE array
implicit none
! Arrays to store the rhs, for time integration
- COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments_rhs ! (ipj,ikr,ikz,updatetlevel)
+ COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_rhs_e ! (ip,ij,ikr,ikz,updatetlevel)
+ COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_rhs_i ! (ip,ij,ikr,ikz,updatetlevel)
! Intermediate steps in rhs of equations
!COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE:: moments_Apl, moments_Bpl, moments_Cpl, moments_Dpl,&
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 0dcbedafd5018ffdc1d211685e7ea89a8af6f652..20e71f3ae116598a7dca5066d352934236c660ec 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -18,6 +18,7 @@ MODULE basic
INTEGER :: iframe1d ! counting the number of times 1d datasets are outputed (for diagnose)
INTEGER :: iframe2d ! counting the number of times 2d datasets are outputed (for diagnose)
INTEGER :: iframe3d ! counting the number of times 3d datasets are outputed (for diagnose)
+ INTEGER :: iframe5d ! counting the number of times 5d datasets are outputed (for diagnose)
! List of logical file units
INTEGER :: lu_in = 90 ! File duplicated from STDIN
@@ -25,7 +26,7 @@ MODULE basic
INTERFACE allocate_array
MODULE PROCEDURE allocate_array_dp1,allocate_array_dp2,allocate_array_dp3,allocate_array_dp4
- MODULE PROCEDURE allocate_array_dc1,allocate_array_dc2,allocate_array_dc3,allocate_array_dc4
+ MODULE PROCEDURE allocate_array_dc1,allocate_array_dc2,allocate_array_dc3,allocate_array_dc4, allocate_array_dc5
MODULE PROCEDURE allocate_array_i1,allocate_array_i2,allocate_array_i3,allocate_array_i4
MODULE PROCEDURE allocate_array_l1,allocate_array_l2,allocate_array_l3,allocate_array_l4
END INTERFACE allocate_array
@@ -67,6 +68,7 @@ CONTAINS
! To allocate arrays of doubles, integers, etc. at run time
SUBROUTINE allocate_array_dp1(a,is1,ie1)
+ IMPLICIT NONE
real(dp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1
ALLOCATE(a(is1:ie1))
@@ -74,6 +76,7 @@ CONTAINS
END SUBROUTINE allocate_array_dp1
SUBROUTINE allocate_array_dp2(a,is1,ie1,is2,ie2)
+ IMPLICIT NONE
real(dp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
ALLOCATE(a(is1:ie1,is2:ie2))
@@ -81,6 +84,7 @@ CONTAINS
END SUBROUTINE allocate_array_dp2
SUBROUTINE allocate_array_dp3(a,is1,ie1,is2,ie2,is3,ie3)
+ IMPLICIT NONE
real(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3))
@@ -88,6 +92,7 @@ CONTAINS
END SUBROUTINE allocate_array_dp3
SUBROUTINE allocate_array_dp4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
+ IMPLICIT NONE
real(dp), DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4))
@@ -95,6 +100,7 @@ CONTAINS
END SUBROUTINE allocate_array_dp4
SUBROUTINE allocate_array_dp5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
+ IMPLICIT NONE
real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5))
@@ -104,6 +110,7 @@ CONTAINS
!========================================
SUBROUTINE allocate_array_dc1(a,is1,ie1)
+ IMPLICIT NONE
DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1
ALLOCATE(a(is1:ie1))
@@ -111,6 +118,7 @@ CONTAINS
END SUBROUTINE allocate_array_dc1
SUBROUTINE allocate_array_dc2(a,is1,ie1,is2,ie2)
+ IMPLICIT NONE
DOUBLE COMPLEX, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
ALLOCATE(a(is1:ie1,is2:ie2))
@@ -118,6 +126,7 @@ CONTAINS
END SUBROUTINE allocate_array_dc2
SUBROUTINE allocate_array_dc3(a,is1,ie1,is2,ie2,is3,ie3)
+ IMPLICIT NONE
DOUBLE COMPLEX, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3))
@@ -125,6 +134,7 @@ CONTAINS
END SUBROUTINE allocate_array_dc3
SUBROUTINE allocate_array_dc4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
+ IMPLICIT NONE
DOUBLE COMPLEX, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4))
@@ -132,7 +142,8 @@ CONTAINS
END SUBROUTINE allocate_array_dc4
SUBROUTINE allocate_array_dc5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
- real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ IMPLICIT NONE
+ DOUBLE COMPLEX, DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5))
a=CMPLX(0.0_dp,0.0_dp)
@@ -141,6 +152,7 @@ CONTAINS
!========================================
SUBROUTINE allocate_array_i1(a,is1,ie1)
+ IMPLICIT NONE
INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1
ALLOCATE(a(is1:ie1))
@@ -148,6 +160,7 @@ CONTAINS
END SUBROUTINE allocate_array_i1
SUBROUTINE allocate_array_i2(a,is1,ie1,is2,ie2)
+ IMPLICIT NONE
INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
ALLOCATE(a(is1:ie1,is2:ie2))
@@ -155,6 +168,7 @@ CONTAINS
END SUBROUTINE allocate_array_i2
SUBROUTINE allocate_array_i3(a,is1,ie1,is2,ie2,is3,ie3)
+ IMPLICIT NONE
INTEGER, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3))
@@ -162,6 +176,7 @@ CONTAINS
END SUBROUTINE allocate_array_i3
SUBROUTINE allocate_array_i4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
+ IMPLICIT NONE
INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4))
@@ -169,6 +184,7 @@ CONTAINS
END SUBROUTINE allocate_array_i4
SUBROUTINE allocate_array_i5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
+ IMPLICIT NONE
real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5))
@@ -178,6 +194,7 @@ CONTAINS
!========================================
SUBROUTINE allocate_array_l1(a,is1,ie1)
+ IMPLICIT NONE
LOGICAL, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1
ALLOCATE(a(is1:ie1))
@@ -185,6 +202,7 @@ CONTAINS
END SUBROUTINE allocate_array_l1
SUBROUTINE allocate_array_l2(a,is1,ie1,is2,ie2)
+ IMPLICIT NONE
LOGICAL, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
ALLOCATE(a(is1:ie1,is2:ie2))
@@ -192,6 +210,7 @@ CONTAINS
END SUBROUTINE allocate_array_l2
SUBROUTINE allocate_array_l3(a,is1,ie1,is2,ie2,is3,ie3)
+ IMPLICIT NONE
LOGICAL, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3))
@@ -199,6 +218,7 @@ CONTAINS
END SUBROUTINE allocate_array_l3
SUBROUTINE allocate_array_l4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
+ IMPLICIT NONE
LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4))
@@ -206,10 +226,22 @@ CONTAINS
END SUBROUTINE allocate_array_l4
SUBROUTINE allocate_array_l5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
+ IMPLICIT NONE
real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5))
a=.false.
END SUBROUTINE allocate_array_l5
+ RECURSIVE FUNCTION Factorial(n) RESULT(Fact)
+ IMPLICIT NONE
+ INTEGER :: Fact
+ INTEGER, INTENT(IN) :: n
+ IF (n == 0) THEN
+ Fact = 1
+ ELSE
+ Fact = n * Factorial(n-1)
+ END IF
+ END FUNCTION Factorial
+
END MODULE basic
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 0a8f93a1768b9977f0966fff74f7a7e8af9778be..7e8c419f543cb19761ada2af1d397aa66b6cd269 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -50,7 +50,7 @@ SUBROUTINE diagnose(kstep)
!CALL creatg(fidres, "/data/var0d", "0d history arrays")
!CALL creatg(fidres, "/data/var1d", "1d profiles")
CALL creatg(fidres, "/data/var2d", "2d profiles")
- CALL creatg(fidres, "/data/var3d", "3d profiles")
+ CALL creatg(fidres, "/data/var5d", "5d profiles")
! Initialize counter of number of saves for each category
@@ -63,15 +63,15 @@ SUBROUTINE diagnose(kstep)
END IF
CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
IF (cstep==0) THEN
- iframe3d=0
+ iframe5d=0
END IF
- CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
+ CALL attach(fidres,"/data/var5d/" , "frames", iframe3d)
! File group
CALL creatg(fidres, "/files", "files")
CALL attach(fidres, "/files", "jobnum", jobnum)
- ! var2d group
+ ! var2d group (electro. pot.)
rank = 0
CALL creatd(fidres, rank, dims, "/data/var2d/time", "Time t*c_s/R")
CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number")
@@ -81,15 +81,22 @@ SUBROUTINE diagnose(kstep)
CALL putarr(fidres, "/data/var2d/phi/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
END IF
- ! var3d group
+ ! var5d group (moments)
rank = 0
- CALL creatd(fidres, rank, dims, "/data/var3d/time", "Time t*c_s/R")
- CALL creatd(fidres, rank, dims, "/data/var3d/cstep", "iteration number")
+ CALL creatd(fidres, rank, dims, "/data/var5d/time", "Time t*c_s/R")
+ CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number")
IF (write_moments) THEN
- CALL creatg(fidres, "/data/var3d/moments", "moments")
- CALL putarr(fidres, "/data/var3d/moments/coordpj", pjarray(ipjs:ipje),"(Jmaxa+1)*p+j+1",ionode=0)
- CALL putarr(fidres, "/data/var3d/moments/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
- CALL putarr(fidres, "/data/var3d/moments/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+ CALL creatg(fidres, "/data/var5d/moments_e", "moments_e")
+ CALL putarr(fidres, "/data/var5d/moments_e/coordp", parray_e(ips_e:ipe_e), "p_e",ionode=0)
+ CALL putarr(fidres, "/data/var5d/moments_e/coordj", jarray_e(ijs_e:ije_e), "j_e",ionode=0)
+ CALL putarr(fidres, "/data/var5d/moments_e/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var5d/moments_e/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+
+ CALL creatg(fidres, "/data/var5d/moments_i", "moments_i")
+ CALL putarr(fidres, "/data/var5d/moments_i/coordp", parray_i(ips_i:ipe_i), "p_i",ionode=0)
+ CALL putarr(fidres, "/data/var5d/moments_i/coordj", jarray_i(ijs_i:ije_i), "j_i",ionode=0)
+ CALL putarr(fidres, "/data/var5d/moments_i/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var5d/moments_i/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
END IF
! Add input namelist variables as attributes of /data/input, defined in srcinfo.h
@@ -145,7 +152,7 @@ SUBROUTINE diagnose(kstep)
!________________________________________________________________________________
! 2. Periodic diagnostics
!
- IF (kstep .GT. 0 .OR. kstep .EQ. 0) THEN
+ IF (kstep .GE. 0) THEN
! 2.1 0d history arrays
IF (nsave_0d .NE. 0) THEN
@@ -165,9 +172,9 @@ SUBROUTINE diagnose(kstep)
END IF
! 2.4 3d profiles
- IF (nsave_3d .NE. 0) THEN
- IF (MOD(cstep, nsave_3d) == 0) THEN
- CALL diagnose_3d
+ IF (nsave_5d .NE. 0) THEN
+ IF (MOD(cstep, nsave_5d) == 0) THEN
+ CALL diagnose_5d
END IF
END IF
@@ -245,45 +252,65 @@ CONTAINS
END SUBROUTINE write_field2d
END SUBROUTINE diagnose_2d
-
-SUBROUTINE diagnose_3d
+
+ SUBROUTINE diagnose_5d
USE basic
USE futils, ONLY: append, getatt, attach, putarrnd
USE fields
+ USE fourier_grid, only: ips_e,ipe_e, ips_i, ipe_i, &
+ ijs_e,ije_e, ijs_i, ije_i
USE time_integration
USE diagnostics_par
use prec_const
IMPLICIT NONE
- CALL append(fidres, "/data/var3d/time", time,ionode=0)
- CALL append(fidres, "/data/var3d/cstep", real(cstep,dp),ionode=0)
- CALL getatt(fidres, "/data/var3d/", "frames",iframe3d)
- iframe3d=iframe3d+1
- CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
+ CALL append(fidres, "/data/var5d/time", time,ionode=0)
+ CALL append(fidres, "/data/var5d/cstep", real(cstep,dp),ionode=0)
+ CALL getatt(fidres, "/data/var5d/", "frames",iframe5d)
+ iframe5d=iframe5d+1
+ CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
IF (write_moments) THEN
- CALL write_field3d(moments(:,:,:,updatetlevel), 'moments')
+ CALL write_field5d_e(moments_e(:,:,:,:,updatetlevel), 'moments_e')
+ CALL write_field5d_i(moments_i(:,:,:,:,updatetlevel), 'moments_i')
END IF
CONTAINS
- SUBROUTINE write_field3d(field, text)
- USE futils, ONLY: attach, putarr
- USE fourier_grid, only: ipjs,ipje, ikrs,ikre, ikzs,ikze
- use prec_const
+ SUBROUTINE write_field5d_e(field, text)
+ USE futils, ONLY: attach, putarr
+ USE fourier_grid, only: ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze
+ USE prec_const
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(ipjs:ipje,ikrs:ikre,ikzs:ikze), INTENT(IN) :: field
+ COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
CHARACTER(LEN=50) :: dset_name
- WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
- CALL putarr(fidres, dset_name, field(ipjs:ipje,ikrs:ikre,ikzs:ikze),ionode=0)
+ WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
+ CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze),ionode=0)
CALL attach(fidres, dset_name, "time", time)
- END SUBROUTINE write_field3d
-
- END SUBROUTINE diagnose_3d
\ No newline at end of file
+ END SUBROUTINE write_field5d_e
+
+ SUBROUTINE write_field5d_i(field, text)
+ USE futils, ONLY: attach, putarr
+ USE fourier_grid, only: ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze
+ USE prec_const
+ IMPLICIT NONE
+
+ COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze), INTENT(IN) :: field
+ CHARACTER(*), INTENT(IN) :: text
+
+ CHARACTER(LEN=50) :: dset_name
+
+ WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
+ CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze),ionode=0)
+
+ CALL attach(fidres, dset_name, "time", time)
+
+ END SUBROUTINE write_field5d_i
+ END SUBROUTINE diagnose_5d
\ No newline at end of file
diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90
index 075bff14fc0611a2b77ec94b1aaa5bd6f5b0fea2..39de0ab7e443533d7fbded74962d94439eca8b89 100644
--- a/src/diagnostics_par_mod.F90
+++ b/src/diagnostics_par_mod.F90
@@ -9,8 +9,7 @@ MODULE diagnostics_par
LOGICAL, PUBLIC, PROTECTED :: write_phi=.TRUE.
LOGICAL, PUBLIC, PROTECTED :: write_doubleprecision=.FALSE.
- INTEGER, PUBLIC, PROTECTED :: nsave_0d , nsave_1d , nsave_2d , nsave_3d
- ! REAL(dp), PUBLIC :: last_timeout_0d, last_timeout_1d, last_timeout_2d, last_timeout_3d
+ INTEGER, PUBLIC, PROTECTED :: nsave_0d , nsave_1d , nsave_2d , nsave_5d
! HDF5 file
@@ -33,7 +32,7 @@ CONTAINS
USE prec_const
IMPLICIT NONE
- NAMELIST /OUTPUT_PAR/ nsave_0d , nsave_1d , nsave_2d , nsave_3d
+ NAMELIST /OUTPUT_PAR/ nsave_0d , nsave_1d , nsave_2d , nsave_5d
NAMELIST /OUTPUT_PAR/ write_moments, write_phi, write_doubleprecision
NAMELIST /OUTPUT_PAR/ resfile0!, rstfile0
@@ -59,7 +58,7 @@ CONTAINS
CALL attach(fidres, TRIM(str), "nsave_0d", nsave_0d)
CALL attach(fidres, TRIM(str), "nsave_1d", nsave_1d)
CALL attach(fidres, TRIM(str), "nsave_2d", nsave_2d)
- CALL attach(fidres, TRIM(str), "nsave_3d", nsave_3d)
+ CALL attach(fidres, TRIM(str), "nsave_5d", nsave_5d)
CALL attach(fidres, TRIM(str), "resfile0", resfile0)
END SUBROUTINE output_par_outputinputs
diff --git a/src/fields_mod.F90 b/src/fields_mod.F90
index 00f8271986bfd0e4e88e975203d98144bed79e33..2e97e6b69f78123f9d974780499b6b22ecb0fad5 100644
--- a/src/fields_mod.F90
+++ b/src/fields_mod.F90
@@ -3,8 +3,9 @@ MODULE fields
use prec_const
implicit none
!------------------MOMENTS Napj------------------
- ! Hermite-Moments: N_a^pj ! dimensions correspond to: ipj, kr, kz, updatetlevel.
- COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments
+ ! Hermite-Moments: N_a^pj ! dimensions correspond to: p, j, kr, kz, updatetlevel.
+ COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_e
+ COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_i
!------------------ELECTROSTATIC POTENTIAL------------------
diff --git a/src/fourier_grid_mod.F90 b/src/fourier_grid_mod.F90
index abf12e7df69d674b5aa02cbfc0582ceed4e78105..c8a40656c5f238e61031832de0e64697ed5ef473 100644
--- a/src/fourier_grid_mod.F90
+++ b/src/fourier_grid_mod.F90
@@ -17,9 +17,9 @@ MODULE fourier_grid
REAL(dp), PUBLIC, PROTECTED :: kzmin = 0._dp ! kz coordinate for left boundary
REAL(dp), PUBLIC, PROTECTED :: kzmax = 1._dp ! kz coordinate for right boundary
- ! Indices of s -> start , e-> end
- INTEGER, PUBLIC, PROTECTED :: ipjs, ipje
- INTEGER, PUBLIC, PROTECTED :: Nmome, Nmomi, Nmomtot
+ ! Indices of s -> start , e-> END
+ INTEGER, PUBLIC, PROTECTED :: ips_e,ipe_e, ijs_e,ije_e
+ INTEGER, PUBLIC, PROTECTED :: ips_i,ipe_i, ijs_i,ije_i
INTEGER, PUBLIC, PROTECTED :: ikrs, ikre, ikzs, ikze
! Toroidal direction
@@ -31,66 +31,67 @@ MODULE fourier_grid
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: krarray
! Grid containing the polynomials degrees
- INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: pjarray
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: parray_e
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: parray_i
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: jarray_e
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: jarray_i
! Public Functions
PUBLIC :: set_krgrid, set_kzgrid, set_pj
PUBLIC :: fourier_grid_readinputs, fourier_grid_outputinputs
- PUBLIC :: bare, bari
- PUBLIC :: rabe, rabi
-contains
+CONTAINS
- subroutine set_krgrid
- use prec_const
- implicit none
- integer :: ikr
- ! Start and end indices of grid
+ SUBROUTINE set_krgrid
+ USE prec_const
+ IMPLICIT NONE
+ INTEGER :: ikr
+ ! Start and END indices of grid
ikrs = 1
ikre = nkr
! Grid spacings, precompute some inverses
- deltakr = (krmax - krmin) / real(nkr,dp)
+ deltakr = (krmax - krmin) / REAL(nkr,dp)
! Discretized kr positions
- allocate(krarray(ikrs:ikre))
+ ALLOCATE(krarray(ikrs:ikre))
DO ikr = ikrs,ikre
- krarray(ikr) = krmin + real(ikr-1,dp) * deltakr
+ krarray(ikr) = krmin + REAL(ikr-1,dp) * deltakr
END DO
- end subroutine set_krgrid
+ END SUBROUTINE set_krgrid
- subroutine set_kzgrid
- use prec_const
- implicit none
- integer :: ikz
- ! Start and end indices of grid
+ SUBROUTINE set_kzgrid
+ USE prec_const
+ IMPLICIT NONE
+ INTEGER :: ikz
+ ! Start and END indices of grid
ikzs = 1
ikze = nkz
! Grid spacings, precompute some inverses
- deltakz = (kzmax - kzmin) / real(nkz,dp)
+ deltakz = (kzmax - kzmin) / REAL(nkz,dp)
! Discretized kz positions
- allocate(kzarray(ikzs:ikze))
+ ALLOCATE(kzarray(ikzs:ikze))
DO ikz = ikzs,ikze
- kzarray(ikz) = kzmin + real(ikz-1,dp) * deltakz
- END DO
- end subroutine set_kzgrid
-
- subroutine set_pj
- use prec_const
- implicit none
- integer :: ipj
- ! number of electrons moments
- Nmome = (Pmaxe + 1)*(Jmaxe + 1)
- ! number of ions moments
- Nmomi = (Pmaxi + 1)*(Jmaxi + 1)
- ! total number of moments
- Nmomtot = Nmome + Nmomi
- ipjs = 1
- ipje = Nmomtot
- ! Polynomials degrees pj = (Jmaxs + 1)*p + j + 1
- allocate(pjarray(ipjs:ipje))
- DO ipj = ipjs,ipje
- pjarray(ipj) = ipj
+ kzarray(ikz) = kzmin + REAL(ikz-1,dp) * deltakz
END DO
- end subroutine set_pj
+ END SUBROUTINE set_kzgrid
+
+ SUBROUTINE set_pj
+ USE prec_const
+ IMPLICIT NONE
+ INTEGER :: ip, ij
+ ips_e = 1; ipe_e = pmaxe + 1
+ ips_i = 1; ipe_i = pmaxi + 1
+ ALLOCATE(parray_e(ips_e:ipe_e))
+ ALLOCATE(parray_i(ips_i:ipe_i))
+ DO ip = ips_e,ipe_e; parray_e(ip) = ip-1; END DO
+ DO ip = ips_i,ipe_i; parray_i(ip) = ip-1; END DO
+
+ ijs_e = 1; ije_e = jmaxe + 1
+ ijs_i = 1; ije_i = jmaxi + 1
+ ALLOCATE(jarray_e(ijs_e:ije_e))
+ ALLOCATE(jarray_i(ijs_i:ije_i))
+ DO ij = ijs_e,ije_e; jarray_e(ij) = ij-1; END DO
+ DO ij = ijs_i,ije_i; jarray_i(ij) = ij-1; END DO
+ END SUBROUTINE set_pj
SUBROUTINE fourier_grid_readinputs
! Read the input parameters
@@ -131,37 +132,4 @@ contains
CALL attach(fidres, TRIM(str), "kzmax", kzmax)
END SUBROUTINE fourier_grid_outputinputs
- !============To handle p,j coefficients efficiently
- FUNCTION bare(p,j) RESULT(idx)
- USE prec_const
- IMPLICIT NONE
- INTEGER, INTENT(in) :: p,j
- INTEGER :: idx
-
- idx = (jmaxe + 1)*p + j + 1
-
- END FUNCTION bare
-
- FUNCTION bari(p,j) RESULT(idx)
- INTEGER, INTENT(in) :: p,j
- INTEGER :: idx
-
- idx = Nmome + (jmaxi + 1)*p + j + 1
-
- END FUNCTION bari
-
- FUNCTION rabe(idx) RESULT(pj)
- INTEGER, INTENT(in) :: idx
- INTEGER, DIMENSION(2) :: pj
- pj(1) = int(FLOOR(real(idx-1) / (jmaxe + 1)))
- pj(2) = idx - 1 - pj(1) * (jmaxe+1)
- END FUNCTION rabe
-
- FUNCTION rabi(idx) RESULT(pj)
- INTEGER, INTENT (in):: idx
- INTEGER, DIMENSION(2) :: pj
- pj(1) = FLOOR(real((idx-Nmome - 1) / (jmaxi + 1)))
- pj(2) = (idx-Nmome) - 1 - pj(1) * (jmaxi+1)
- END FUNCTION rabi
-
END MODULE fourier_grid
diff --git a/src/inital.F90 b/src/inital.F90
index 2699fc6c57e7829878c60f05bc48aefabb80315c..49e8f9a8b4c7d3a8fcc4d4174607fe25e416f2e6 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -24,7 +24,7 @@ SUBROUTINE init_profiles
use prec_const
IMPLICIT NONE
- INTEGER :: ipj, ikr, ikz
+ INTEGER :: ip,ij, ikr,ikz
INTEGER, DIMENSION(12) :: iseedarr
REAL(dp) :: noise
@@ -36,11 +36,26 @@ SUBROUTINE init_profiles
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
- DO ipj=ipjs,ipje
+
+ DO ip=ips_e,ipe_e
+ DO ij=ijs_e,ije_e
+ CALL RANDOM_NUMBER(noise)
+ moments_e( ip,ij, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) ! Re
+ CALL RANDOM_NUMBER(noise)
+ moments_e( ip,ij, ikr,ikz, :) = imagu * (initback_moments + initnoise_moments*(noise-0.5_dp)) ! Im
+ END DO
+ END DO
+
+ DO ip=ips_i,ipe_i
+ DO ij=ijs_i,ije_i
+ CALL RANDOM_NUMBER(noise)
+ moments_i( ip,ij, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) ! Re
CALL RANDOM_NUMBER(noise)
- moments( ipj, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+ moments_i( ip,ij, ikr,ikz, :) = imagu * (initback_moments + initnoise_moments*(noise-0.5_dp)) ! Im
END DO
END DO
+
+ END DO
END DO
CALL poisson ! To set phi
diff --git a/src/memory.F90 b/src/memory.F90
index 035697f0528d56ddfc4bb2a0e947b4a31ec5ac67..af5933d0be418d4640d608491677f041311ef3fd 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -10,22 +10,13 @@ SUBROUTINE memory
use prec_const
IMPLICIT NONE
-
-
! Moments and moments rhs
- CALL allocate_array( moments, ipjs,ipje, ikrs,ikre, ikzs,ikze, 1,ntimelevel)
- CALL allocate_array( moments_rhs, ipjs,ipje, ikrs,ikre, ikzs,ikze, 1,ntimelevel)
+ CALL allocate_array( moments_e, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze, 1,ntimelevel )
+ CALL allocate_array( moments_i, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze, 1,ntimelevel )
+ CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze, 1,ntimelevel )
+ CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze, 1,ntimelevel )
! Electrostatic potential
- CALL allocate_array( phi, ikrs,ikre, ikzs,ikze)
+ CALL allocate_array(phi, ikrs,ikre, ikzs,ikze)
- ! Intermediate steps in rhs of equations
- !CALL allocate_array( moments_Apl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
- !CALL allocate_array( moments_Bpl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
- !CALL allocate_array( moments_Cpl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
- !CALL allocate_array( moments_Dpl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
- !CALL allocate_array( moments_Epl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
- !CALL allocate_array( moments_Fpl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
- !CALL allocate_array( moments_Gpl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
- !CALL allocate_array( moments_Hpl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
END SUBROUTINE memory
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index f2178b30be2a0521866fc8a2666b99bb58325e67..58171853c3a908687868072971e0da974ce7ccbc 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -10,99 +10,254 @@ SUBROUTINE moments_eq_rhs
use prec_const
IMPLICIT NONE
- INTEGER :: ip,ij, ipj, ipp2j, ipm2j, ipjp1, ipjm1, ikr,ikz
- INTEGER :: kronp0, kronp2
- INTEGER, DIMENSION(2) :: tmppj
+ INTEGER :: ip,ij, ikr,ikz ! loops indices
REAL(dp) :: ip_dp, ij_dp
- REAL(dp) :: kz, taueqeetaBi, tauiqietaBi, tauaqaetaBi
- COMPLEX(dp) :: Napj, Napp2j, Napm2j, Napjp1, Napjm1
-
- taueqeetaBi = tau_e/real(q_e)*eta_B * imagu ! Factor tau_s/qu_s*eta_B*i
- tauiqietaBi = tau_i/real(q_i)*eta_B * imagu
-
- Napp2j = 0; Napm2j = 0; Napjp1 = 0; Napjm1 = 0; ! higher mixing moments terms
- kronp0 = 0; kronp2 = 0; ! Useful kronecker for phi terms
-
- ipp2j = 1; ipm2j = 1; ipjp1 = 1; ! Mixing moment indices are set to one initially
- ipjm1 = 1; ! in order to not overpass the moment array
-
- do ipj = ipjs, ipje
-
- if (ipj .le. Nmome + 1) then ! electrons moments
-
- tmppj = rabe(ipj) ! Compute p,j electrons moments degree
- ip = tmppj(1); ij = tmppj(2);
- if (ip+2 .le. pmaxe+1) then
- ipp2j = bare(ip+2,ij)
- Napp2j = moments(ipp2j,ikr,ikz,updatetlevel)
- endif
- if (ip-2 .ge. 1) then
- ipm2j = bare(ip-2,ij)
- Napm2j = moments(ipm2j,ikr,ikz,updatetlevel)
- endif
- if (ij+1 .le. jmaxe+1) then
- ipjp1 = bare(ip,ij+1)
- Napjp1 = moments(ipjp1,ikr,ikz,updatetlevel)
- endif
- if (ij-1 .ge. 1) then
- ipjm1 = bare(ip,ij-1)
- Napjm1 = moments(ipjm1,ikr,ikz,updatetlevel)
- endif
- tauaqaetaBi = taueqeetaBi !species dependant factor
-
- else ! ions moments
-
- tmppj = rabi(ipj) ! Compute p,j ions moments degree
- ip = tmppj(1); ij = tmppj(2);
- if (ip+2 .le. pmaxi+1) then
- ipp2j = bari(ip+2,ij)
- Napp2j = moments(ipp2j,ikr,ikz,updatetlevel)
- endif
- if (ip-2 .ge. 1) then
- ipm2j = bari(ip-2,ij)
- Napm2j = moments(ipm2j,ikr,ikz,updatetlevel)
- endif
- if (ij+1 .le. jmaxi+1) then
- ipjp1 = bari(ip,ij+1)
- Napjp1 = moments(ipjp1,ikr,ikz,updatetlevel)
- endif
- if (ij-1 .ge. 1) then
- ipjm1 = bari(ip,ij-1)
- Napjm1 = moments(ipjm1,ikr,ikz,updatetlevel)
- endif
- tauaqaetaBi = tauiqietaBi
- endif
-
- if (ip .eq. 0) then
- kronp0 = 1
- endif
- if (ip .eq. 2) then
- kronp2 = 1
- endif
-
- Napj = moments(ipj,ikr,ikz,updatetlevel)
-
- ip_dp = real(ip,dp)
- ij_dp = real(ij,dp)
-
- do ikr = ikrs,ikre
- do ikz = ikzs,ikze
- kz = kzarray(ikz)
- ! N_a^{pj} term
- Napj = 2*(ip_dp + ij_dp + 1._dp) * Napj
- ! N_a^{p+2,j} term
- Napp2j = SQRT(ip_dp + 1._dp) * SQRT(ip_dp + 2._dp) * Napp2j
- ! N_a^{p-2,j} term
- Napm2j = SQRT(ip_dp) * SQRT(ip_dp - 1._dp) * Napm2j
- ! N_a^{p,j+1} term
- Napjp1 = -(ij_dp + 1._dp) * Napjp1
- ! N_a^{p,j-1} term
- Napjm1 = -ij_dp * Napjm1
-
- moments_rhs(ipj,ikr,ikz,updatetlevel) = &
- tauaqaetaBi * kz * (Napj + Napp2j + Napm2j + Napjp1 + Napjm1)
- end do
- end do
- end do
+ COMPLEX(dp) :: kr, kz, kperp2
+ COMPLEX(dp) :: taue_qe_etaBi, taui_qi_etaBi
+ COMPLEX(dp) :: kernelj, kerneljp1, kerneljm1, b_e2, b_i2 ! Kernel functions and variable
+ COMPLEX(dp) :: factj, xb_e2, xb_i2 ! Auxiliary variables
+ COMPLEX(dp) :: xNapj, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! factors depending on the pj loop
+ COMPLEX(dp) :: xphij, xphijp1, xphijm1, xColl
+ COMPLEX(dp) :: TNapj, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi, TColl ! terms of the rhs
+
+ !!!!!!!!! Electrons moments RHS !!!!!!!!!
+ Tphi = 0 ! electrostatic potential term
+ taue_qe_etaBi = tau_e/q_e * eta_B * imagu ! one-time computable factor
+ xb_e2 = sigma_e**2 * tau_e/2. ! species dependant factor of the Kernel, squared
+
+ ploope : DO ip = ips_e, ipe_e
+ ip_dp = real(ip-1,dp) ! real index is one minus the loop index (0 to pmax)
+
+ ! x N_e^{p+2,j}
+ IF (ip+2 .LE. pmaxe + 1) THEN
+ xNapp2j = taue_qe_etaBi * SQRT(ip_dp + 1._dp) * SQRT(ip_dp + 2._dp)
+ ELSE
+ xNapp2j = 0.
+ ENDIF
+
+ ! x N_e^{p-2,j}
+ IF (ip-2 .GE. 1) THEN
+ xNapm2j = taue_qe_etaBi * SQRT(ip_dp) * SQRT(ip_dp - 1.)
+ ELSE
+ xNapm2j = 0.
+ ENDIF
+
+ jloope : DO ij = ijs_e, ije_e
+ ij_dp = real(ij-1,dp) ! real index is one minus the loop index (0 to jmax)
+
+ ! x N_e^{p,j+1}
+ IF (ij+1 .LE. jmaxe+1) THEN
+ xNapjp1 = -taue_qe_etaBi * (ij_dp + 1.)
+ ELSE
+ xNapjp1 = 0.
+ ENDIF
+
+ ! x N_e^{p,j-1}
+ IF (ij-1 .GE. 1) THEN
+ xNapjm1 = -taue_qe_etaBi * ij_dp
+ ELSE
+ xNapjm1 = 0.
+ ENDIF
+
+ ! x N_e^{pj}
+ xNapj = taue_qe_etaBi * 2.*(ip_dp + ij_dp + 1.)
+
+ ! first part of collision operator (Lenard-Bernstein)
+ xColl = -nu*(ip_dp + 2.*ij_dp)
+
+ ! x phi
+ IF (ip .eq. 1) THEN !(krokecker delta_p^0)
+ xphij = imagu * (eta_n + 2.*ij_dp*eta_T + 2.*eta_B*(ij_dp+1.) )
+ xphijp1 = -imagu * (eta_B + eta_T)*(ij_dp+1.)
+ xphijm1 = -imagu * (eta_B + eta_T)* ij_dp
+ factj = real(Factorial(ij-1),dp)
+ ELSE IF (ip .eq. 3) THEN !(krokecker delta_p^2)
+ xphij = imagu * (SQRT2*eta_B + eta_T/SQRT2)
+ factj = real(Factorial(ij-1),dp)
+ ELSE
+ xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
+ factj = 1
+ ENDIF
+
+ ! Loop on kspace
+ krloope : DO ikr = ikrs,ikre
+ kzloope : DO ikz = ikzs,ikze
+ kr = krarray(ikr) ! Poloidal wavevector
+ kz = kzarray(ikz) ! Toroidal wavevector
+ kperp2 = kr**2 + kz**2 ! perpendicular wavevector
+ b_e2 = kperp2 * xb_e2 ! Bessel argument
+
+ !! Compute moments and mixing terms
+ ! term propto N_e^{p,j}
+ TNapj = moments_e(ip,ij,ikr,ikz,updatetlevel) * xNapj
+ ! term propto N_e^{p+2,j}
+ IF (xNapp2j .NE. (0.,0.)) THEN
+ TNapp2j = moments_e(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j
+ ELSE
+ TNapp2j = 0.
+ ENDIF
+ ! term propto N_e^{p-2,j}
+ IF (xNapm2j .NE. (0.,0.)) THEN
+ TNapm2j = moments_e(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j
+ ELSE
+ TNapm2j = 0.
+ ENDIF
+ ! xterm propto N_e^{p,j+1}
+ IF (xNapjp1 .NE. (0.,0.)) THEN
+ TNapjp1 = moments_e(ip,ij+1,ikr,ikz,updatetlevel) * xNapp2j
+ ELSE
+ TNapjp1 = 0.
+ ENDIF
+ ! term propto N_e^{p,j-1}
+ IF (xNapjm1 .NE. (0.,0.)) THEN
+ TNapjm1 = moments_e(ip,ij-1,ikr,ikz,updatetlevel) * xNapp2j
+ ELSE
+ TNapjm1 = 0.
+ ENDIF
+
+ ! Collision term completed (Lenard-Bernstein)
+ IF (xNapp2j .NE. (0.,0.)) THEN
+ TColl = (xColl - nu * b_e2/4.) * moments_e(ip+2,ij,ikr,ikz,updatetlevel)
+ ELSE
+ TColl = 0.
+ ENDIF
+
+ !! Electrical potential term
+ Tphi = 0
+ IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! 0 otherwise (krokecker delta_p^0)
+ kernelj = b_e2**(ij-1) * exp(-b_e2)/factj
+ kerneljp1 = kernelj * b_e2 /(ij_dp + 1.)
+ kerneljm1 = kernelj * ij_dp / b_e2
+ Tphi = (xphij * Kernelj + xphijp1 * Kerneljp1 + xphijm1 * Kerneljm1)* kz * phi(ikr,ikz)
+ ENDIF
+
+ ! Sum of all precomputed terms
+ moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
+ kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1) &
+ + Tphi + TColl
+
+ END DO kzloope
+ END DO krloope
+
+ END DO jloope
+ END DO ploope
+
+
+ !!!!!!!!! Ions moments RHS !!!!!!!!!
+ Tphi = 0 ! electrostatic potential term
+ taui_qi_etaBi = tau_i/real(q_i)*eta_B * imagu ! one-time computable factor
+ xb_i2 = sigma_i**2 * tau_i/2.0 ! species dependant factor of the Kernel, squared
+
+ ploopi : DO ip = ips_i, ipe_i
+ ip_dp = real(ip-1,dp)
+
+ ! x N_i^{p+2,j}
+ IF (ip+2 .LE. pmaxi + 1) THEN
+ xNapp2j = taui_qi_etaBi * SQRT(ip_dp + 1.) * SQRT(ip_dp + 2.)
+ ELSE
+ xNapp2j = 0.
+ ENDIF
+
+ ! x N_i^{p-2,j}
+ IF (ip-2 .GE. 1) THEN
+ xNapm2j = taui_qi_etaBi * SQRT(ip_dp) * SQRT(ip_dp - 1.)
+ ELSE
+ xNapm2j = 0.
+ ENDIF
+
+ jloopi : DO ij = ijs_i, ije_i
+ ij_dp = real(ij-1,dp)
+
+ ! x N_i^{p,j+1}
+ IF (ij+1 .LE. jmaxi+1) THEN
+ xNapjp1 = -taui_qi_etaBi * (ij_dp + 1.)
+ ELSE
+ xNapjp1 = 0.
+ ENDIF
+
+ ! x N_i^{p,j-1}
+ IF (ij-1 .GE. 1) THEN
+ xNapjm1 = -taui_qi_etaBi * ij_dp
+ ELSE
+ xNapjm1 = 0.
+ ENDIF
+
+ ! x N_i^{pj}
+ xNapj = taui_qi_etaBi * 2.*(ip_dp + ij_dp + 1.)
+
+ ! first part of collision operator (Lenard-Bernstein)
+ xColl = -nu*(ip_dp + 2.0*ij_dp)
+
+ ! x phi
+ IF (ip .EQ. 1) THEN !(krokecker delta_p^0)
+ xphij = imagu * (eta_n + 2.*ij_dp*eta_T + 2.*eta_B*(ij_dp+1.) )
+ xphijp1 = -imagu * (eta_B + eta_T)*(ij_dp+1.)
+ xphijm1 = -imagu * (eta_B + eta_T)* ij_dp
+ factj = REAL(Factorial(ij-1),dp)
+ ELSE IF (ip .EQ. 3) THEN !(krokecker delta_p^2)
+ xphij = imagu * (SQRT2*eta_B + eta_T/SQRT2)
+ factj = REAL(Factorial(ij-1),dp)
+ ELSE
+ xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
+ ENDIF
+ ! Loop on kspace
+ krloopi : DO ikr = ikrs,ikre
+ kzloopi : DO ikz = ikzs,ikze
+ kr = krarray(ikr) ! Poloidal wavevector
+ kz = kzarray(ikz) ! Toroidal wavevector
+ kperp2 = kr**2 + kz**2 ! perpendicular wavevector
+ b_i2 = kperp2 * xb_i2 ! Bessel argument
+
+ !! Compute moments and mixing terms
+ ! term propto N_i^{p,j}
+ TNapj = moments_i(ip,ij,ikr,ikz,updatetlevel) * xNapj
+
+ IF (xNapp2j .NE. (0.,0.)) THEN ! term propto N_i^{p+2,j}
+ TNapp2j = moments_i(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j
+ ELSE
+ TNapp2j = 0.
+ ENDIF
+ IF (xNapm2j .NE. (0.,0.)) THEN ! term propto N_i^{p-2,j}
+ TNapm2j = moments_i(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j
+ ELSE
+ TNapm2j = 0.
+ ENDIF
+ IF (xNapjp1 .NE. (0.,0.)) THEN ! xterm propto N_a^{p,j+1}
+ TNapjp1 = moments_i(ip,ij+1,ikr,ikz,updatetlevel) * xNapp2j
+ ELSE
+ TNapjp1 = 0.
+ ENDIF
+ IF (xNapjm1 .NE. (0.,0.)) THEN ! term propto N_a^{p,j-1}
+ TNapjm1 = moments_i(ip,ij-1,ikr,ikz,updatetlevel) * xNapp2j
+ ELSE
+ TNapjm1 = 0.
+ ENDIF
+
+ ! Collision term completed (Lenard-Bernstein)
+ IF (xNapp2j .NE. (0.,0.)) THEN
+ TColl = (xColl - nu * b_i2/4.) * moments_i(ip+2,ij,ikr,ikz,updatetlevel)
+ ELSE
+ TColl = 0.
+ ENDIF
+
+ !! Electrical potential term
+ Tphi = 0
+ IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! 0 otherwise (krokecker delta_p^0, delta_p^2)
+ kernelj = b_i2**(ij-1) * exp(-b_i2)/factj
+ kerneljp1 = kernelj * b_i2 /(ij_dp + 1._dp)
+ kerneljm1 = kernelj * ij_dp / b_i2
+ Tphi = (xphij * Kernelj + xphijp1 * Kerneljp1 + xphijm1 * Kerneljm1)* kz * phi(ikr,ikz)
+ ENDIF
+
+ ! Sum of all precomputed terms
+ moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
+ kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1) &
+ + Tphi + TColl
+
+ END DO kzloopi
+ END DO krloopi
+
+ END DO jloopi
+ END DO ploopi
END SUBROUTINE moments_eq_rhs
diff --git a/src/poisson.F90 b/src/poisson.F90
index cae10fdc46933758757cc4f4b8c0ff99cb3ec47d..f7ab0c5886495c8e3af252eab4507ea51204cf75 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -1,5 +1,4 @@
-
-subroutine poisson
+SUBROUTINE poisson
! Solve poisson equation to get phi
USE time_integration, ONLY: updatetlevel
@@ -12,58 +11,73 @@ subroutine poisson
IMPLICIT NONE
INTEGER :: ikr,ikz, ini,ine, i0j
- REAL(dp) :: kr, kz
- REAL(dp) :: k1_, k2_ ! sub kernel factor for recursive build
- REAL(dp) :: x_e, x_i, alphaD
+ REAL(dp) :: ini_dp, ine_dp
+ COMPLEX(dp) :: kr, kz, kperp2
+ COMPLEX(dp) :: xkm_, xk2_ ! sub kernel factor for recursive build
+ COMPLEX(dp) :: b_e2, b_i2, alphaD
+ COMPLEX(dp) :: sigmaetaue2, sigmaitaui2 ! To avoid redondant computation
COMPLEX(dp) :: gammaD, gammaphi
COMPLEX(dp) :: sum_kernel2_e, sum_kernel2_i
COMPLEX(dp) :: sum_kernel_mom_e, sum_kernel_mom_i
+ sigmaetaue2 = sigma_e**2 * tau_e/2.
+ sigmaitaui2 = sigma_i**2 * tau_i/2.
+
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
- kr = krarray(ikr)
- kz = kzarray(ikz)
-
- ! Compute electrons sum(Kernel**2) and sum(Kernel * Ne0j)
- x_e = sqrt(kr**2 + kz**2) * sigma_e * sqrt(2.0*tau_e)/2.0
- sum_kernel2_e = 0
- sum_kernel_mom_e = 0
- k1_ = moments(1,ikr,ikz,updatetlevel)
- k2_ = 1.0
- if (jmaxe .ge. 0) then
- DO ine=1,jmaxe
- i0j = bare(0,ine)
- k1_ = k1_ * x_e**2/ine
- k2_ = k2_ * x_e**4/ine**2
- sum_kernel_mom_e = sum_kernel_mom_e + k1_ * moments(i0j,ikr,ikz,updatetlevel)
- sum_kernel2_e = sum_kernel2_e + k2_
+ kr = krarray(ikr)
+ kz = kzarray(ikz)
+ kperp2 = kr**2 + kz**2
+
+ ! Compute electrons sum(Kernel * Ne0n) (skm) and sum(Kernel**2) (sk2)
+ b_e2 = kperp2 * sigmaetaue2 ! non dim kernel argument (kperp2 sigma_a sqrt(2 tau_a)/2)
+
+ xkm_ = 1. ! Initialization for n = 0
+ xk2_ = 1.
+
+ sum_kernel_mom_e = moments_e(1,1,ikr,ikz,updatetlevel)
+ sum_kernel2_e = xk2_
+
+ if (jmaxe .GT. 0) then
+ DO ine=2,jmaxe+1
+ ine_dp = REAL(ine-1,dp)
+
+ xkm_ = xkm_ * b_e2/ine_dp
+ xk2_ = xk2_ *(b_e2/ine_dp)**2
+
+ sum_kernel_mom_e = sum_kernel_mom_e + xkm_ * moments_e(1,ine,ikr,ikz,updatetlevel)
+ sum_kernel2_e = sum_kernel2_e + xk2_
END DO
endif
- sum_kernel2_e = sum_kernel2_e * exp(-x_e**2)
- sum_kernel_mom_e = sum_kernel_mom_e * exp(-x_e**2)**2
-
- ! Compute ions sum(Kernel**2) and sum(Kernel * Ne0j)
- x_i = sqrt(kr**2 + kz**2) * sigma_i * sqrt(2.0*tau_i)/2.0
- sum_kernel2_i = 0
- sum_kernel_mom_i = 0
- k1_ = moments(Nmomi + 1, ikr, ikz, updatetlevel)
- k2_ = 1.0
- if (jmaxi .ge. 0) then
- DO ini=1,jmaxe
- i0j = bari(0,ini)
- k1_ = k1_ * x_i**2/ini
- k2_ = k2_ * x_i**4/ini**2
- sum_kernel_mom_i = sum_kernel_mom_i + k1_ * moments(Nmome + i0j,ikr,ikz,updatetlevel)
- sum_kernel2_i = sum_kernel2_i + k2_
+ sum_kernel2_e = sum_kernel2_e * exp(-b_e2)
+ sum_kernel_mom_e = sum_kernel_mom_e * exp(-2.*b_e2)
+
+ ! Compute ions sum(Kernel * Ni0n) (skm) and sum(Kernel**2) (sk2)
+ b_i2 = kperp2 * sigmaitaui2
+
+ xkm_ = 1.
+ xk2_ = 1.
+
+ sum_kernel_mom_i = moments_i(1, 1, ikr, ikz, updatetlevel)
+ sum_kernel2_i = xk2_
+ if (jmaxi .GT. 0) then
+ DO ini=2,jmaxi + 1
+ ini_dp = REAL(ini-1,dp) ! Real index (0 to jmax)
+
+ xkm_ = xkm_ * b_i2/ini_dp
+ xk2_ = xk2_ *(b_i2/ini_dp)**2
+
+ sum_kernel_mom_i = sum_kernel_mom_i + xkm_ * moments_i(1,ini,ikr,ikz,updatetlevel)
+ sum_kernel2_i = sum_kernel2_i + xk2_
END DO
endif
- sum_kernel2_i = sum_kernel2_i * exp(-x_i**2)
- sum_kernel_mom_i = sum_kernel_mom_i * exp(-x_i**2)**2
+ sum_kernel2_i = sum_kernel2_i * exp(-b_i2)
+ sum_kernel_mom_i = sum_kernel_mom_i * exp(-2.*b_i2)
! Assembling the poisson equation
- alphaD = sqrt(kr**2 + kz**2) * lambdaD**2
- gammaD = alphaD + q_e**2/tau_e * sum_kernel2_e &
- + q_i**2/tau_i * sum_kernel2_i
+ alphaD = kperp2 * lambdaD**2.
+ gammaD = alphaD + (q_e**2)/tau_e * sum_kernel2_e &
+ + (q_i**2)/tau_i * sum_kernel2_i
gammaphi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i
@@ -72,4 +86,4 @@ subroutine poisson
END DO
END DO
-END SUBROUTINE poisson
+END SUBROUTINE poisson
\ No newline at end of file
diff --git a/src/stepon.F90 b/src/stepon.F90
index db13a317a8cc16167e0252d7c2dbe7c7a9149f24..bbd75a2bd101e0c65c1613d185e15edd8bddd354 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -3,8 +3,8 @@ SUBROUTINE stepon
USE basic
USE time_integration
- USE fields, ONLY: moments, phi
- USE array , ONLY: moments_rhs
+ USE fields, ONLY: moments_e, moments_i, phi
+ USE array , ONLY: moments_rhs_e, moments_rhs_i
USE fourier_grid
USE advance_field_routine, ONLY: advance_time_level, advance_field
USE model
@@ -13,21 +13,34 @@ SUBROUTINE stepon
use prec_const
IMPLICIT NONE
- INTEGER :: num_step, ipj
+ INTEGER :: num_step, ip,ij
DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
! Compute right hand side
+ !WRITE (*,*) 'Compute right hand side .. nstep = ', num_step
CALL moments_eq_rhs
+ !WRITE (*,*) 'Advance time level .. nstep = ', num_step
CALL advance_time_level ! Advance from updatetlevel to updatetlevel+1
- do ipj=ipjs,ipje
- CALL advance_field(moments(ipj,:,:,:),moments_rhs(ipj,:,:,:))
- enddo
+ !WRITE (*,*) 'Advance field .. nstep = ', num_step
+
+ DO ip=ips_e,ipe_e
+ DO ij=ijs_e,ije_e
+ CALL advance_field(moments_e(ip,ij,:,:,:),moments_rhs_e(ip,ij,:,:,:))
+ ENDDO
+ ENDDO
+
+ DO ip=ips_i,ipe_i
+ DO ij=ijs_i,ije_i
+ CALL advance_field(moments_i(ip,ij,:,:,:),moments_rhs_i(ip,ij,:,:,:))
+ ENDDO
+ ENDDO
! Solving Poisson equation
+ !WRITE (*,*) 'Solving Poisson equation .. nstep = ', num_step
CALL poisson
CALL checkfield_all()
@@ -38,9 +51,17 @@ SUBROUTINE stepon
SUBROUTINE checkfield_all ! Check all the fields for inf or nan
IF(.NOT.nlend) THEN
nlend=nlend .or. checkfield(phi,' phi')
- do ipj=ipjs,ipje
- nlend=nlend .or. checkfield(moments(ipj,:,:,updatetlevel),' moments')
- enddo
+ DO ip=ips_e,ipe_e
+ DO ij=ijs_e,ije_e
+ nlend=nlend .or. checkfield(moments_e(ip,ij,:,:,updatetlevel),' moments_e')
+ ENDDO
+ ENDDO
+
+ DO ip=ips_i,ipe_i
+ DO ij=ijs_i,ije_i
+ nlend=nlend .or. checkfield(moments_i(ip,ij,:,:,updatetlevel),' moments_i')
+ ENDDO
+ ENDDO
ENDIF
END SUBROUTINE checkfield_all
diff --git a/wk/fort.90 b/wk/fort.90
index 46dbdca373ac1d12fa11f30d82161380a8a7a001..bc01967580ff4f8cbef78e9ad3cc58eb934ad302 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -3,51 +3,51 @@ A Skeleton for MPI Time-Dependent program
S. Brunner, T.M. Tran CRPP/EPFL
-
&BASIC
- nrun=3
- dt=1.0D-3
+ nrun=10000
+ dt=0.01
tmax=100 ! time normalized to 1/omega_pe
/
&GRID
- pmaxe=10 ! number of Hermite moments (including Ne,Te,Ue)
- jmaxe=5 ! number of Hermite moments (including Ne,Te,Ue)
- pmaxi=10 ! number of Hermite moments (including Ne,Te,Ue)
- jmaxi=5 ! number of Hermite moments (including Ne,Te,Ue)
- nkr=1
- krmin=0.
- krmax=0. ! Normalized to ?
- nkz=1
- kzmin=1.
- kzmax=1. ! Normalized to ?
+ pmaxe = 15 ! number of Hermite moments
+ jmaxe = 4 ! number of Hermite moments
+ pmaxi = 15 ! number of Hermite moments
+ jmaxi = 4 ! number of Hermite moments
+ nkr = 1
+ krmin = 0.
+ krmax = 0. ! Normalized to ?
+ nkz = 1
+ kzmin = 0.1
+ kzmax = 0.1 ! Normalized to ?
/
&OUTPUT_PAR
- nsave_0d=0
- nsave_1d=0
- nsave_2d=1
- nsave_3d=1
- write_moments=.true.
- write_phi=.true.
- write_doubleprecision=.true.
- resfile0='results'
+ nsave_0d = 0
+ nsave_1d = 0
+ nsave_2d = 100
+ nsave_5d = 100
+ write_moments = .true.
+ write_phi = .true.
+ write_doubleprecision = .true.
+ resfile0 = 'results'
/
&MODEL_PAR
! Collisionality
- nu = 0.01 ! Normalized collision frequency normalized to plasma frequency
+ nu = 0.001 ! Normalized collision frequency normalized to plasma frequency
tau_e = 1.0 ! T_e/T_e
tau_i = 1.0 ! T_i/T_e temperature ratio
sigma_e = 0.0233380 ! sqrt(m_e/m_i) mass ratio
sigma_i = 1.0 ! sqrt(m_i/m_i)
q_e =-1.0 ! Electrons charge
q_i = 1.0 ! Ions charge
- eta_n = 1.0 ! L_perp / L_n Density gradient
+ eta_n = 2.0 ! L_perp / L_n Density gradient
eta_T = 0.0 ! L_perp / L_T Temperature gradient
- eta_B = 0.5 ! L_perp / L_B Magnetic gradient and curvature
+ eta_B = 1.0 ! L_perp / L_B Magnetic gradient and curvature
lambdaD = 0.0 ! Debye length
/
&INITIAL_CON
! Background value
- initback_moments=1. ! x 1e-3
+ initback_moments=0.01 ! x 1e-3
! Noise amplitude
- initnoise_moments=0.1
+ initnoise_moments=0.
iseed=42
/
&TIME_INTEGRATION_PAR