diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
index d30b3a8c5299564a0a9ff20b8b29a8e514a4d3f1..1ef419a8c7dd0440da7ea0dbd845c1bef52c2c0e 100644
--- a/src/calculus_mod.F90
+++ b/src/calculus_mod.F90
@@ -20,25 +20,25 @@ MODULE calculus
CONTAINS
-SUBROUTINE grad_z(target,local_Nz,Ngz,inv_deltaz,f,ddzf)
+SUBROUTINE grad_z(target,local_nz,Ngz,inv_deltaz,f,ddzf)
implicit none
! Compute the periodic boundary condition 4 points centered finite differences
! formula among staggered grid or not.
! not staggered : the derivative results must be on the same grid as the field
! staggered : the derivative is computed from a grid to the other
- INTEGER, INTENT(IN) :: target, local_Nz, Ngz
+ INTEGER, INTENT(IN) :: target, local_nz, Ngz
REAL(dp), INTENT(IN) :: inv_deltaz
- COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: f
- COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: ddzf
+ COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: f
+ COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddzf
INTEGER :: iz
IF(Ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
SELECT CASE(TARGET)
CASE(1)
- CALL grad_z_o2e(local_Nz,Ngz,inv_deltaz,f,ddzf)
+ CALL grad_z_o2e(local_nz,Ngz,inv_deltaz,f,ddzf)
CASE(2)
- CALL grad_z_e2o(local_Nz,Ngz,inv_deltaz,f,ddzf)
+ CALL grad_z_e2o(local_nz,Ngz,inv_deltaz,f,ddzf)
CASE DEFAULT ! No staggered grid -> usual centered finite differences
- DO iz = 1,local_Nz
+ DO iz = 1,local_nz
ddzf(iz) = dz_usu(-2)*f(iz ) + dz_usu(-1)*f(iz+1) &
+dz_usu( 0)*f(iz+2) &
+dz_usu( 1)*f(iz+3) + dz_usu( 2)*f(iz+4)
@@ -49,30 +49,30 @@ SUBROUTINE grad_z(target,local_Nz,Ngz,inv_deltaz,f,ddzf)
ddzf = 0._dp
ENDIF
CONTAINS
- SUBROUTINE grad_z_o2e(local_Nz,Ngz,inv_deltaz,fo,ddzfe) ! Paruta 2018 eq (27)
+ SUBROUTINE grad_z_o2e(local_nz,Ngz,inv_deltaz,fo,ddzfe) ! Paruta 2018 eq (27)
! gives the gradient of a field from the odd grid to the even one
implicit none
- INTEGER, INTENT(IN) :: local_Nz, Ngz
+ INTEGER, INTENT(IN) :: local_nz, Ngz
REAL(dp), INTENT(IN) :: inv_deltaz
- COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: fo
- COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: ddzfe !
+ COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: fo
+ COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddzfe !
INTEGER :: iz
- DO iz = 1,local_Nz
+ DO iz = 1,local_nz
ddzfe(iz) = dz_o2e(-2)*fo(iz ) + dz_o2e(-1)*fo(iz+1) &
+dz_o2e( 0)*fo(iz+2) + dz_o2e( 1)*fo(iz+3)
ENDDO
ddzfe(:) = ddzfe(:) * inv_deltaz
END SUBROUTINE grad_z_o2e
- SUBROUTINE grad_z_e2o(local_Nz,Ngz,inv_deltaz,fe,ddzfo) ! n2v for Paruta 2018 eq (28)
+ SUBROUTINE grad_z_e2o(local_nz,Ngz,inv_deltaz,fe,ddzfo) ! n2v for Paruta 2018 eq (28)
! gives the gradient of a field from the even grid to the odd one
implicit none
- INTEGER, INTENT(IN) :: local_Nz, Ngz
+ INTEGER, INTENT(IN) :: local_nz, Ngz
REAL(dp), INTENT(IN) :: inv_deltaz
- COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: fe
- COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: ddzfo
+ COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: fe
+ COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddzfo
INTEGER :: iz
- DO iz = 1,local_Nz
+ DO iz = 1,local_nz
ddzfo(iz) = dz_e2o(-1)*fe(iz+1) + dz_e2o(0)*fe(iz+2) &
+dz_e2o( 1)*fe(iz+3) + dz_e2o(2)*fe(iz+4)
ENDDO
@@ -80,16 +80,16 @@ CONTAINS
END SUBROUTINE grad_z_e2o
END SUBROUTINE grad_z
-SUBROUTINE grad_z2(local_Nz,Ngz,inv_deltaz,f,ddz2f)
+SUBROUTINE grad_z2(local_nz,Ngz,inv_deltaz,f,ddz2f)
! Compute the second order fourth derivative for periodic boundary condition
implicit none
- INTEGER, INTENT(IN) :: local_Nz, Ngz
+ INTEGER, INTENT(IN) :: local_nz, Ngz
REAL(dp), INTENT(IN) :: inv_deltaz
- COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: f
- COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: ddz2f
+ COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: f
+ COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddz2f
INTEGER :: iz
IF(Ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
- DO iz = 1,local_Nz
+ DO iz = 1,local_nz
ddz2f(iz) = dz2_usu(-2)*f(iz ) + dz2_usu(-1)*f(iz+1) &
+dz2_usu( 0)*f(iz+2)&
+dz2_usu( 1)*f(iz+3) + dz2_usu( 2)*f(iz+4)
@@ -101,16 +101,16 @@ SUBROUTINE grad_z2(local_Nz,Ngz,inv_deltaz,f,ddz2f)
END SUBROUTINE grad_z2
-SUBROUTINE grad_z4(local_Nz,Ngz,inv_deltaz,f,ddz4f)
+SUBROUTINE grad_z4(local_nz,Ngz,inv_deltaz,f,ddz4f)
! Compute the second order fourth derivative for periodic boundary condition
implicit none
- INTEGER, INTENT(IN) :: local_Nz, Ngz
+ INTEGER, INTENT(IN) :: local_nz, Ngz
REAL(dp), INTENT(IN) :: inv_deltaz
- COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: f
- COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: ddz4f
+ COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: f
+ COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddz4f
INTEGER :: iz
IF(Ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
- DO iz = 1,local_Nz
+ DO iz = 1,local_nz
ddz4f(iz) = dz4_usu(-2)*f(iz ) + dz4_usu(-1)*f(iz+1) &
+dz4_usu( 0)*f(iz+2)&
+dz4_usu( 1)*f(iz+3) + dz4_usu( 2)*f(iz+4)
@@ -122,47 +122,47 @@ SUBROUTINE grad_z4(local_Nz,Ngz,inv_deltaz,f,ddz4f)
END SUBROUTINE grad_z4
-SUBROUTINE interp_z(target,local_Nz,Ngz,f_in,f_out)
+SUBROUTINE interp_z(target,local_nz,Ngz,f_in,f_out)
! Function meant to interpolate one field defined on a even/odd z into
! the other odd/even z grid.
! If Staggered Grid flag (SG) is false, returns identity
implicit none
- INTEGER, INTENT(IN) :: local_Nz, Ngz
+ INTEGER, INTENT(IN) :: local_nz, Ngz
INTEGER, intent(in) :: target ! target grid : 0 for even grid, 1 for odd
- COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: f_in
- COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: f_out
+ COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: f_in
+ COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: f_out
SELECT CASE(TARGET)
CASE(1) ! output on even grid
- CALL interp_o2e_z(local_Nz,Ngz,f_in,f_out)
+ CALL interp_o2e_z(local_nz,Ngz,f_in,f_out)
CASE(2) ! output on odd grid
- CALL interp_e2o_z(local_Nz,Ngz,f_in,f_out)
+ CALL interp_e2o_z(local_nz,Ngz,f_in,f_out)
CASE DEFAULT ! No staggered grid -> usual centered finite differences
f_out = f_in
END SELECT
CONTAINS
- SUBROUTINE interp_o2e_z(local_Nz, Ngz,fo,fe)
+ SUBROUTINE interp_o2e_z(local_nz, Ngz,fo,fe)
! gives the value of a field from the odd grid to the even one
implicit none
- INTEGER, INTENT(IN) :: local_Nz, Ngz
- COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: fo
- COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: fe
+ INTEGER, INTENT(IN) :: local_nz, Ngz
+ COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: fo
+ COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: fe
INTEGER :: iz
! 4th order interp
- DO iz = 1,local_Nz
+ DO iz = 1,local_nz
fe(iz) = iz_o2e(-2)*fo(iz ) + iz_o2e(-1)*fo(iz+1) &
+ iz_o2e( 0)*fo(iz+2) + iz_o2e( 1)*fo(iz+3)
ENDDO
END SUBROUTINE interp_o2e_z
- SUBROUTINE interp_e2o_z(local_Nz, Ngz,fe,fo)
+ SUBROUTINE interp_e2o_z(local_nz, Ngz,fe,fo)
! gives the value of a field from the even grid to the odd one
implicit none
- INTEGER, INTENT(IN) :: local_Nz, Ngz
- COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: fe
- COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: fo
+ INTEGER, INTENT(IN) :: local_nz, Ngz
+ COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: fe
+ COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: fo
INTEGER :: iz
! 4th order interp
- DO iz = 1,local_Nz
+ DO iz = 1,local_nz
fo(iz) = iz_e2o(-1)*fe(iz+1) + iz_e2o( 0)*fe(iz+2) &
+ iz_e2o( 1)*fe(iz+3) + iz_e2o( 2)*fe(iz+4)
ENDDO
@@ -170,23 +170,23 @@ CONTAINS
END SUBROUTINE interp_z
-SUBROUTINE simpson_rule_z(local_Nz,dz,f,intf)
+SUBROUTINE simpson_rule_z(local_nz,dz,f,intf)
! integrate f(z) over z using the simpon's rule. Assume periodic boundary conditions (f(ize+1) = f(izs))
!from molix BJ Frei
USE prec_const, ONLY: dp, onethird
USE parallel, ONLY: num_procs_z, rank_z, comm_z, manual_0D_bcast
USE mpi
implicit none
- INTEGER, INTENT(IN) :: local_Nz
+ INTEGER, INTENT(IN) :: local_nz
REAL(dp),INTENT(IN) :: dz
- complex(dp),dimension(local_Nz), intent(in) :: f
+ complex(dp),dimension(local_nz), intent(in) :: f
COMPLEX(dp), intent(out) :: intf
COMPLEX(dp) :: buffer, local_int
INTEGER :: root, i_, iz, ierr
! Buil local sum using the weights of composite Simpson's rule
local_int = 0._dp
- DO iz = 1,local_Nz
+ DO iz = 1,local_nz
IF(MODULO(iz,2) .EQ. 1) THEN ! odd iz
local_int = local_int + 2._dp*onethird*dz*f(iz)
ELSE ! even iz
diff --git a/src/control.F90 b/src/control.F90
index 944daf4cfc3bf7758c77e2022e3ae89189a22328..74e41c9d3e38ce9ddec854de88d0e59ca1dc6615 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -61,6 +61,9 @@ SUBROUTINE control
! 2. Main loop
DO
CALL cpu_time(t0_step) ! Measuring time
+
+ CALL tesend
+ IF( nlend ) EXIT ! exit do loop
CALL increase_step
CALL increase_cstep
@@ -68,8 +71,6 @@ SUBROUTINE control
CALL increase_time
- CALL tesend
- IF( nlend ) EXIT ! exit do loop
CALL diagnose(step)
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index a98d5981c5c8218f0274ea89e940c7322e08f2af..f1a7d8ee2ee525143d51bafb6b0047bed7484543 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -1,9 +1,9 @@
SUBROUTINE diagnose(kstep)
! Diagnostics, writing simulation state to disk
- USE basic, ONLY: t0_diag,t1_diag,tc_diag, lu_in, finish, start, cstep, dt, time, tmax, display_h_min_s
+ USE basic, ONLY: t0_diag,t1_diag,tc_diag, lu_in, finish, start, cstep, dt, time, tmax, display_h_min_s
USE diagnostics_par, ONLY: input_fname
- USE processing, ONLY: pflux_x, hflux_x
- USE parallel, ONLY: my_id
+ USE processing, ONLY: pflux_x, hflux_x
+ USE parallel, ONLY: my_id
IMPLICIT NONE
INTEGER, INTENT(in) :: kstep
CALL cpu_time(t0_diag) ! Measuring time
@@ -24,7 +24,7 @@ SUBROUTINE diagnose(kstep)
!! Specific diagnostic calls
CALL diagnose_full(kstep)
! Terminal info
- IF ((kstep .GT. 0) .AND. (MOD(cstep, INT(1.0/dt)) == 0) .AND. (my_id .EQ. 0)) THEN
+ IF ((kstep .GE. 0) .AND. (MOD(cstep, INT(1.0/dt)) == 0) .AND. (my_id .EQ. 0)) THEN
WRITE(*,"(A,F6.0,A1,F6.0,A8,G10.2,A8,G10.2,A)")'|t/tmax = ', time,"/",tmax,'| Gxi = ',pflux_x(1),'| Qxi = ',hflux_x(1),'|'
ENDIF
CALL cpu_time(t1_diag); tc_diag = tc_diag + (t1_diag - t0_diag)
@@ -155,7 +155,7 @@ SUBROUTINE diagnose_full(kstep)
CALL putarrnd(fidres, "/data/metric/Jacobian", Jacobian(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(1:local_nky,1:local_nkx,1+ngz/2:local_nz+ngz/2,:), (/1, 1, 3/))
- CALL putarrnd(fidres, "/data/metric/kernel", kernel(1,1+ngj/2:local_nj+ngj/2,1:local_nky,1:local_nkx,1+ngz/2:local_nz+ngz/2,:), (/1, 1, 2, 4/))
+ CALL putarrnd(fidres, "/data/metric/kernel", kernel(1,1+ngj/2:local_nj+ngj/2,1:local_nky,1:local_nkx,1+ngz/2:local_nz+ngz/2,1), (/1, 1, 2, 4/))
! var0d group (gyro transport)
IF (nsave_0d .GT. 0) THEN
CALL creatg(fidres, "/data/var0d", "0d profiles")
@@ -225,13 +225,10 @@ SUBROUTINE diagnose_full(kstep)
CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
END IF
ENDIF
-
-
!_____________________________________________________________________________
! 2. Periodic diagnostics
!
IF (kstep .GE. 0) THEN
-
! 2.1 0d history arrays
IF (nsave_0d .GT. 0) THEN
IF ( MOD(cstep, nsave_0d) == 0 ) THEN
@@ -248,15 +245,13 @@ SUBROUTINE diagnose_full(kstep)
ENDIF
ENDIF
! 2.4 5d profiles
- IF (nsave_5d .GT. 0 .AND. cstep .GT. 0) THEN
+ IF (nsave_5d .GT. 0) THEN
IF (MOD(cstep, nsave_5d) == 0) THEN
CALL diagnose_5d
END IF
END IF
-
!_____________________________________________________________________________
! 3. Final diagnostics
-
ELSEIF (kstep .EQ. -1) THEN
CALL attach(fidres, "/data/input","cpu_time",finish-start)
@@ -446,8 +441,8 @@ SUBROUTINE diagnose_5d
USE fields, ONLY: moments
USE grid, ONLY:total_np, total_nj, total_nky, total_nkx, total_nz, &
local_np, local_nj, local_nky, local_nkx, local_nz, &
- ngp, ngj, ngz
- USE time_integration, ONLY: updatetlevel
+ ngp, ngj, ngz, total_na
+ USE time_integration, ONLY: updatetlevel, ntimelevel
USE diagnostics_par
USE prec_const, ONLY: dp
IMPLICIT NONE
@@ -470,19 +465,19 @@ SUBROUTINE diagnose_5d
USE parallel, ONLY: gather_pjxyz, num_procs
USE prec_const, ONLY: dp
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN) :: field
+ COMPLEX(dp), DIMENSION(total_na,local_np+ngp,local_nj+ngj,local_nky,local_nkx,local_nz+ngz,ntimelevel), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
- COMPLEX(dp), DIMENSION(local_np,local_nj,local_nky,local_nkx,local_nz) :: field_sub
- COMPLEX(dp), DIMENSION(total_np,total_nj,total_nky,total_nkx,total_nz) :: field_full
+ COMPLEX(dp), DIMENSION(total_na,local_np,local_nj,local_nky,local_nkx,local_nz) :: field_sub
+ COMPLEX(dp), DIMENSION(total_na,total_np,total_nj,total_nky,total_nkx,total_nz) :: field_full
CHARACTER(LEN=50) :: dset_name
- field_sub = field(1,(1+ngp/2):(local_np+ngp/2),(1+ngj/2):(local_nj+ngj/2),&
+ field_sub = field(1:total_na,(1+ngp/2):(local_np+ngp/2),(1+ngj/2):(local_nj+ngj/2),&
1:local_nky,1:local_nkx, (1+ngz/2):(local_nz+ngz/2),updatetlevel)
field_full = 0;
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
IF (num_procs .EQ. 1) THEN
CALL putarr(fidres, dset_name, field_sub, ionode=0)
ELSEIF(GATHERV_OUTPUT) THEN ! output using one node (gatherv)
- CALL gather_pjxyz(field_sub,field_full,local_np,total_np,total_nj,local_nky,total_nky,total_nkx,local_nz,total_nz)
+ CALL gather_pjxyz(field_sub,field_full,total_na,local_np,total_np,total_nj,local_nky,total_nky,total_nkx,local_nz,total_nz)
CALL putarr(fidres, dset_name, field_full, ionode=0)
ELSE
CALL putarrnd(fidres, dset_name, field_sub, (/1,3,5/))
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 0af30081e067d1d586b51d9e9ee165a067060b26..cb59aa62757314df65375268a04171081851c32d 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -100,7 +100,7 @@ CONTAINS
! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient
implicit none
REAL(dp) :: kx,ky
- COMPLEX(dp), DIMENSION(local_Nz+Ngz) :: integrant
+ COMPLEX(dp), DIMENSION(local_nz) :: integrant
real(dp) :: G1,G2,G3,Cx,Cy
INTEGER :: eo,iz,iky,ikx
@@ -138,7 +138,7 @@ CONTAINS
CALL set_kparray(gxx,gxy,gyy,hatB)
DO eo = 1,Nzgrid
! Curvature operator (Frei et al. 2022 eq 2.15)
- DO iz = 1,local_Nz+Ngz
+ DO iz = 1,local_nz+Ngz
G1 = gxx(iz,eo)*gyy(iz,eo)-gxy(iz,eo)*gxy(iz,eo)
G2 = gxx(iz,eo)*gyz(iz,eo)-gxy(iz,eo)*gxz(iz,eo)
G3 = gxy(iz,eo)*gyz(iz,eo)-gyy(iz,eo)*gxz(iz,eo)
@@ -166,7 +166,7 @@ CONTAINS
CALL set_ikx_zBC_map
!
! Compute the inverse z integrated Jacobian (useful for flux averaging)
- integrant = Jacobian(:,1) ! Convert into complex array
+ integrant = Jacobian(1+ngz/2:local_nz+ngz/2,1) ! Convert into complex array
CALL simpson_rule_z(local_nz,deltaz,integrant,iInt_Jacobian)
iInt_Jacobian = 1._dp/iInt_Jacobian ! reverse it
END SUBROUTINE eval_magnetic_geometry
@@ -175,7 +175,7 @@ CONTAINS
!
SUBROUTINE eval_salpha_geometry
- USE grid, ONLY : local_Nz,Ngz,zarray,Nzgrid
+ USE grid, ONLY : local_nz,Ngz,zarray,Nzgrid
! evaluate s-alpha geometry model
implicit none
REAL(dp) :: z
@@ -183,7 +183,7 @@ CONTAINS
alpha_MHD = 0._dp
DO eo = 1,Nzgrid
- DO iz = 1,local_Nz+Ngz
+ DO iz = 1,local_nz+Ngz
z = zarray(iz,eo)
! metric
@@ -228,14 +228,14 @@ CONTAINS
!
SUBROUTINE eval_zpinch_geometry
- USE grid, ONLY : local_Nz,Ngz,zarray,Nzgrid
+ USE grid, ONLY : local_nz,Ngz,zarray,Nzgrid
implicit none
REAL(dp) :: z
INTEGER :: iz, eo
alpha_MHD = 0._dp
DO eo = 1,Nzgrid
- DO iz = 1,local_Nz+Ngz
+ DO iz = 1,local_nz+Ngz
z = zarray(iz,eo)
! metric
@@ -278,7 +278,7 @@ CONTAINS
!--------------------------------------------------------------------------------
! NOT TESTED
subroutine eval_1D_geometry
- USE grid, ONLY : local_Nz,Ngz,zarray, Nzgrid
+ USE grid, ONLY : local_nz,Ngz,zarray, Nzgrid
! evaluate 1D perp geometry model
implicit none
REAL(dp) :: z
@@ -313,7 +313,7 @@ CONTAINS
USE grid, ONLY: local_nky,Nkx, contains_zmin,contains_zmax, Nexc
USE prec_const, ONLY: imagu, pi
IMPLICIT NONE
- ! REAL :: shift
+ ! REAL(dp) :: shift
INTEGER :: ikx,iky
ALLOCATE(ikx_zBC_L(local_nky,Nkx))
ALLOCATE(ikx_zBC_R(local_nky,Nkx))
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 75a2a5240ee0315d28b33e2d8f3ae9fb1f25fdec..5b5794fdcc2003d700ad8a7fef0acd2ba4b24944 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -76,7 +76,7 @@ MODULE grid
REAL(dp), PUBLIC, PROTECTED :: local_kxmin, local_kxmax
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: local_zmin, local_zmax
! local z weights for computing simpson rule
- INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: zweights_SR
+ REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: zweights_SR
! Numerical diffusion scaling
REAL(dp), PUBLIC, PROTECTED :: diff_p_coeff, diff_j_coeff
REAL(dp), PUBLIC, PROTECTED :: diff_kx_coeff, diff_ky_coeff, diff_dz_coeff
@@ -381,7 +381,7 @@ CONTAINS
CHARACTER(len=*), INTENT(IN) ::LINEARITY
INTEGER, INTENT(IN) :: N_HD
INTEGER :: ikx, ikxo
- REAL :: Lx_adapted
+ REAL(dp):: Lx_adapted
IF(shear .GT. 0) THEN
IF(my_id.EQ.0) write(*,*) 'Magnetic shear detected: set up sheared kx grid..'
! mininal size of box in x to respect dkx = 2pi shear dky
@@ -424,7 +424,7 @@ CONTAINS
local_kxmax = 0._dp
DO ikx = ikxs,ikxe
ikxo = ikx - local_nkx_offset
- kxarray(ikxo) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1)/real(Nkx)))
+ kxarray(ikxo) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1,dp)/real(Nkx,dp)))
if (ikx .EQ. Nx/2+1) kxarray(ikxo) = -kxarray(ikxo)
! Finding kx=0
IF (kxarray(ikxo) .EQ. 0) THEN
@@ -441,7 +441,7 @@ CONTAINS
! Build the full grids on process 0 to diagnose it without comm
! kx
DO ikx = 1,Nkx
- kxarray_full(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1)/real(Nkx)))
+ kxarray_full(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1,dp)/real(Nkx,dp)))
IF (ikx .EQ. Nx/2+1) kxarray_full(ikx) = -kxarray_full(ikx)
END DO
ELSE ! Odd number of kx (-2 -1 0 1 2)
@@ -503,8 +503,8 @@ CONTAINS
USE prec_const
USE parallel, ONLY: num_procs_z, rank_z
IMPLICIT NONE
- REAL :: grid_shift, Lz, zmax, zmin
- INTEGER :: istart, iend, in, Npol, iz, ig, eo, izo
+ REAL(dp):: grid_shift, Lz, zmax, zmin
+ INTEGER :: istart, iend, in, Npol, iz, ig, eo
total_nz = Nz
! Length of the flux tube (in ballooning angle)
Lz = 2_dp*pi*Npol
@@ -562,10 +562,9 @@ CONTAINS
! Local z array
ALLOCATE(zarray(local_nz+Ngz,Nzgrid))
!! interior point loop
- DO iz = izs,ize
- izo = iz - local_nz_offset
+ DO iz = 1,total_nz
DO eo = 1,Nzgrid
- zarray(izo,eo) = zarray_full(iz) + (eo-1)*grid_shift
+ zarray(iz+ngz/2,eo) = zarray_full(iz) + REAL(eo-1,dp)*grid_shift
ENDDO
ENDDO
ALLOCATE(local_zmax(Nzgrid),local_zmin(Nzgrid))
@@ -573,16 +572,17 @@ CONTAINS
! Find local extrema
local_zmax(eo) = zarray(local_nz+ngz/2,eo)
local_zmin(eo) = zarray(1+ngz/2,eo)
+ print*, zarray
! Fill the ghosts
- DO ig = 1,ngj/2
- zarray(ig,eo) = local_zmin(eo)-(ngz/2+(ig-1))*deltaz
- zarray(local_nz+ig,eo) = local_zmax(eo)+ig*deltaz
+ DO ig = 1,ngz/2
+ zarray(ig,eo) = local_zmin(eo)-REAL(ngz/2-(ig-1),dp)*deltaz
+ zarray(local_nz+ngz/2+ig,eo) = local_zmax(eo)+REAL(ig,dp)*deltaz
ENDDO
! Set up the flags to know if the process contains the tip and/or the tail
! of the z domain (important for z-boundary condition)
- IF(abs(local_zmin(eo) - (zmin+(eo-1)*grid_shift)) .LT. EPSILON(zmin)) &
+ IF(abs(local_zmin(eo) - (zmin+REAL(eo-1,dp)*grid_shift)) .LT. EPSILON(zmin)) &
contains_zmin = .TRUE.
- IF(abs(local_zmax(eo) - (zmax+(eo-1)*grid_shift)) .LT. EPSILON(zmax)) &
+ IF(abs(local_zmax(eo) - (zmax+REAL(eo-1,dp)*grid_shift)) .LT. EPSILON(zmax)) &
contains_zmax = .TRUE.
ENDDO
IF(mod(Nz,2) .NE. 0 ) THEN
@@ -596,7 +596,7 @@ CONTAINS
REAL(dp) :: kx, ky
CALL allocate_array( kparray, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz, 1,2)
DO eo = 1,Nzgrid
- DO iz = 1,local_Nz+Ngz
+ DO iz = 1,local_nz+Ngz
DO iky = 1,local_nky
ky = kyarray(iky)
DO ikx = 1,local_nkx
diff --git a/src/inital.F90 b/src/inital.F90
index f3bbe1b1d3ed01945a3fd79691e7a8f01dbd0f06..61d7e5f731009d28813ee240a4192922ed5ba536 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -181,7 +181,7 @@ SUBROUTINE init_gyrodens
DO iky=1,local_nky
DO iz=1+ngz/2,local_nz+ngz/2
CALL RANDOM_NUMBER(noise)
- IF ( (parray(ip) .EQ. 1) .AND. (jarray(ij) .EQ. 1) ) THEN
+ IF ( (parray(ip) .EQ. 0) .AND. (jarray(ij) .EQ. 0) ) THEN
moments(ia,ip,ij,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
ELSE
moments(ia,ip,ij,iky,ikx,iz,:) = 0._dp
@@ -196,8 +196,6 @@ SUBROUTINE init_gyrodens
ENDIF
END DO
END DO
- print*, SUM(moments)
- stop
! Putting to zero modes that are not in the 2/3 Orszag rule
IF (LINEARITY .NE. 'linear') THEN
DO ikx=1,total_nkx
@@ -249,7 +247,7 @@ SUBROUTINE init_phi
DO ikx=2,total_nkx/2
phi(iky0,ikx,iz) = phi(iky0,total_nkx+2-ikx,iz)
ENDDO
- phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz)) !origin must be real
+ phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz),dp) !origin must be real
END DO
ENDIF
!**** ensure no previous moments initialization
@@ -309,7 +307,7 @@ SUBROUTINE init_phi_ppj
DO ikx=2,total_nkx/2
phi(iky0,ikx,iz) = phi(iky0,total_nkx+2-ikx,iz)
ENDDO
- phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz)) !origin must be real
+ phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz),dp) !origin must be real
END DO
ENDIF
!**** ensure no previous moments initialization
diff --git a/src/memory.F90 b/src/memory.F90
index 97cf719164b45c33f8931c68c790ac16689408a4..3ee7f91cf71f3311cdfd0ef0869630f3d4fcc93f 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -4,7 +4,7 @@ SUBROUTINE memory
USE array
USE basic, ONLY: allocate_array
USE fields
- USE grid, ONLY: local_Na, local_Np,Ngp ,local_Nj,Ngj, local_nky, local_nkx,local_Nz,Ngz, jmax, pmax
+ USE grid, ONLY: local_Na, local_Np,Ngp ,local_Nj,Ngj, local_nky, local_nkx,local_nz,Ngz, jmax, pmax
USE collision
USE time_integration, ONLY: ntimelevel
USE prec_const
@@ -12,30 +12,30 @@ SUBROUTINE memory
IMPLICIT NONE
! Electrostatic potential
- CALL allocate_array( phi, 1,local_nky, 1,local_nkx, 1,local_Nz+Ngz)
- CALL allocate_array( psi, 1,local_nky, 1,local_nkx, 1,local_Nz+Ngz)
- CALL allocate_array(inv_poisson_op, 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( inv_ampere_op, 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( inv_pol_ion, 1,local_nky, 1,local_nkx, 1,local_Nz)
- ! CALL allocate_array(HF_phi_correction_operator, 1,local_nky, 1,local_nkx, 1,local_Nz)
+ CALL allocate_array( phi, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz)
+ CALL allocate_array( psi, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz)
+ CALL allocate_array(inv_poisson_op, 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( inv_ampere_op, 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( inv_pol_ion, 1,local_nky, 1,local_nkx, 1,local_nz)
+ ! CALL allocate_array(HF_phi_correction_operator, 1,local_nky, 1,local_nkx, 1,local_nz)
!Moments related arrays
- CALL allocate_array( dens, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( upar, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( uper, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( Tpar, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( Tper, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( temp, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( Kernel, 1,local_Na, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_Nz+Ngz, 1,2)
- CALL allocate_array( moments, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_Nz+Ngz, 1,ntimelevel )
- CALL allocate_array( Napjz, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_Nz)
- CALL allocate_array( moments_rhs, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nky, 1,local_nkx, 1,local_Nz, 1,ntimelevel )
- CALL allocate_array( nadiab_moments, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_Nz+Ngz)
- CALL allocate_array( ddz_napj, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( ddzND_Napj, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( interp_napj, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( Capj, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( Sapj, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nky, 1,local_nkx, 1,local_Nz)
+ CALL allocate_array( dens, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( upar, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( uper, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( Tpar, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( Tper, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( temp, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( Kernel, 1,local_Na, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz, 1,2)
+ CALL allocate_array( moments, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz, 1,ntimelevel )
+ CALL allocate_array( Napjz, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nz)
+ CALL allocate_array( moments_rhs, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nky, 1,local_nkx, 1,local_nz, 1,ntimelevel )
+ CALL allocate_array( nadiab_moments, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz)
+ CALL allocate_array( ddz_napj, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( ddzND_Napj, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( interp_napj, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( Capj, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( Sapj, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( xnapj, 1,local_Na, 1,local_Np, 1,local_Nj)
CALL allocate_array( xnapp2j, 1,local_Na, 1,local_Np)
CALL allocate_array( xnapp1j, 1,local_Na, 1,local_Np)
@@ -67,9 +67,9 @@ SUBROUTINE memory
!___________________ 2x5D ARRAYS __________________________
!! Collision matrices
IF (GK_CO) THEN !GK collision matrices (one for each kperp)
- CALL allocate_array( Cab_F, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( Cab_T, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,local_nky, 1,local_nkx, 1,local_Nz)
- CALL allocate_array( Caa, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,local_nky, 1,local_nkx, 1,local_Nz)
+ CALL allocate_array( Cab_F, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( Cab_T, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,local_nky, 1,local_nkx, 1,local_nz)
+ CALL allocate_array( Caa, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,local_nky, 1,local_nkx, 1,local_nz)
ELSE !DK collision matrix (same for every k)
CALL allocate_array( Cab_F, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,1, 1,1, 1,1)
CALL allocate_array( Cab_T, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,1, 1,1, 1,1)
diff --git a/src/miller_mod.F90 b/src/miller_mod.F90
index 688561e3077cd597358c7f43dbc63dede7902923..3ecce924a623f6e7bd0fb46f28bf9757e9fb0a04 100644
--- a/src/miller_mod.F90
+++ b/src/miller_mod.F90
@@ -6,7 +6,7 @@ MODULE miller
USE basic
USE parallel, ONLY: my_id, num_procs_z, nbr_U, nbr_D, comm0
! use coordinates,only: gcoor, get_dzprimedz
- USE grid, ONLY: local_Nky, local_Nkx, local_Nz, Ngz, kyarray, kxarray, zarray, Nz, local_nz_offset
+ USE grid, ONLY: local_Nky, local_Nkx, local_nz, Ngz, kyarray, kxarray, zarray, Nz, local_nz_offset
! use discretization
USE lagrange_interpolation
! use par_in, only: beta, sign_Ip_CW, sign_Bt_CW, Npol
@@ -59,12 +59,12 @@ CONTAINS
real(dp), INTENT(INOUT) :: edge_opt ! alpha mhd
real(dp), INTENT(INOUT) :: dpdx_pm_geom ! amplitude mag. eq. pressure grad.
real(dp), INTENT(INOUT) :: C_y, C_xy
- real(dp), dimension(1:local_Nz+Ngz,1:2), INTENT(INOUT) :: &
+ real(dp), dimension(1:local_nz+Ngz,1:2), INTENT(INOUT) :: &
gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,&
dBdx_,dBdy_,Bfield_,jacobian_,&
dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_, &
gradz_coeff_
- real(dp), dimension(1:local_Nky,1:local_Nkx,1:local_Nz+Ngz,1:2), INTENT(INOUT) :: Ckxky_
+ real(dp), dimension(1:local_Nky,1:local_Nkx,1:local_nz+Ngz,1:2), INTENT(INOUT) :: Ckxky_
INTEGER :: iz, ikx, iky, eo
! No parameter in gyacomo yet
real(dp) :: sign_Ip_CW=1 ! current sign (only normal current)
@@ -477,7 +477,7 @@ CONTAINS
call lag3interp(dxdZ_s,chi_s,np_s,dxdZ_out,chi_out,Nz)
! Fill the interior of the geom arrays with the results
do eo=1,2
- DO iz = 1,local_Nz
+ DO iz = 1,local_nz
gxx_(iz+Ngz/2,eo) = gxx_out(iz-local_nz_offset)
gyy_(iz+Ngz/2,eo) = gyy_out(iz-local_nz_offset)
gxz_(iz+Ngz/2,eo) = gxz_out(iz-local_nz_offset)
@@ -511,7 +511,7 @@ CONTAINS
CALL update_ghosts_z(dxdZ_(:,eo))
! Curvature operator (Frei et al. 2022 eq 2.15)
- DO iz = 1,local_Nz+Ngz
+ DO iz = 1,local_nz+Ngz
G1 = gxy_(iz,eo)*gxy_(iz,eo)-gxx_(iz,eo)*gyy_(iz,eo)
G2 = gxy_(iz,eo)*gxz_(iz,eo)-gxx_(iz,eo)*gyz_(iz,eo)
G3 = gyy_(iz,eo)*gxz_(iz,eo)-gxy_(iz,eo)*gyz_(iz,eo)
@@ -536,10 +536,10 @@ CONTAINS
SUBROUTINE update_ghosts_z(fz_)
IMPLICIT NONE
! INTEGER, INTENT(IN) :: nztot_
- REAL(dp), DIMENSION(1:local_Nz+Ngz), INTENT(INOUT) :: fz_
+ REAL(dp), DIMENSION(1:local_nz+Ngz), INTENT(INOUT) :: fz_
REAL(dp), DIMENSION(-2:2) :: buff
INTEGER :: status(MPI_STATUS_SIZE), count, last, first
- last = local_Nz+Ngz/2
+ last = local_nz+Ngz/2
first= 1 + Ngz/2
IF(Nz .GT. 1) THEN
IF (num_procs_z .GT. 1) THEN
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 189ac0a9d41bb328753275524243cb1c9d8d9852..ff0949bb43ea6ffe6b84878ca2d5312b6c3ee973 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -72,7 +72,8 @@ END SUBROUTINE build_dv4Hp_table
SUBROUTINE evaluate_kernels
USE basic
USE array, ONLY : kernel!, HF_phi_correction_operator
- USE grid, ONLY : local_Na, local_Nj,Ngj, local_nkx, local_nky, local_nz, Ngz, jarray, kparray
+ USE grid, ONLY : local_Na, local_Nj,Ngj, local_nkx, local_nky, local_nz, Ngz, jarray, kparray,&
+ nzgrid
USE species, ONLY : sigma2_tau_o2
USE prec_const, ONLY: dp
IMPLICIT NONE
@@ -80,26 +81,26 @@ SUBROUTINE evaluate_kernels
REAL(dp) :: j_dp, y_, factj
DO ia = 1,local_Na
- DO eo = 1,2
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO iz = 1,local_nz + Ngz
- DO ij = 1, local_Nj + Ngj
- j_int = jarray(ij)
- j_dp = REAL(j_int,dp)
- y_ = sigma2_tau_o2(ia) * kparray(iky,ikx,iz,eo)**2
- IF(j_int .LT. 0) THEN
- kernel(ia,ij,iky,ikx,iz,eo) = 0._dp
- ELSE
- factj = GAMMA(j_dp+1._dp)
- kernel(ia,ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/factj
- ENDIF
+ DO eo = 1,nzgrid
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO iz = 1,local_nz + Ngz
+ DO ij = 1, local_nj + Ngj
+ j_int = jarray(ij)
+ j_dp = REAL(j_int,dp)
+ y_ = sigma2_tau_o2(ia) * kparray(iky,ikx,iz,eo)**2
+ IF(j_int .LT. 0) THEN
+ kernel(ia,ij,iky,ikx,iz,eo) = 0._dp
+ ELSE
+ factj = GAMMA(j_dp+1._dp)
+ kernel(ia,ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/factj
+ ENDIF
+ ENDDO
+ ! IF (ijs .EQ. 1) & ! if ijs is 1, the ghost kernel has negative index
+ ! kernel(ia,ijgs,iky,ikx,iz,eo) = 0._dp
+ ENDDO
+ ENDDO
ENDDO
- ! IF (ijs .EQ. 1) & ! if ijs is 1, the ghost kernel has negative index
- ! kernel(ia,ijgs,iky,ikx,iz,eo) = 0._dp
- ENDDO
- ENDDO
- ENDDO
ENDDO
! !! Correction term for the evaluation of the heat flux
! HF_phi_correction_operator(:,:,:) = &
diff --git a/src/parallel_mod.F90 b/src/parallel_mod.F90
index d90918ef172d99b91088875679fae6d2de349e4a..2cf09efabfc82a79b3186cd0781f11c9eedc0b02 100644
--- a/src/parallel_mod.F90
+++ b/src/parallel_mod.F90
@@ -276,39 +276,39 @@ CONTAINS
!!!!! Gather a field in kinetic + spatial coordinates on rank 0 !!!!!
!!!!! Gather a field in spatial coordinates on rank 0 !!!!!
- SUBROUTINE gather_pjxyz(field_loc,field_tot,np_loc,np_tot,nj_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot)
+ SUBROUTINE gather_pjxyz(field_loc,field_tot,na_tot,np_loc,np_tot,nj_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: np_loc,np_tot,nj_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot
- COMPLEX(dp), DIMENSION(np_loc,nj_tot,nky_loc,nkx_tot,nz_loc), INTENT(IN) :: field_loc
- COMPLEX(dp), DIMENSION(np_tot,nj_tot,nky_tot,nkx_tot,nz_tot), INTENT(OUT) :: field_tot
- COMPLEX(dp), DIMENSION(np_tot,nky_tot,nz_loc) :: buffer_pt_yt_zl ! full p, full y, local z
- COMPLEX(dp), DIMENSION(np_tot,nky_tot,nz_tot) :: buffer_pt_yt_zt ! full p, full y, full z
- COMPLEX(dp), DIMENSION(np_tot,nky_loc) :: buffer_pt_yl_zc ! full p, local y, constant z
- COMPLEX(dp), DIMENSION(np_tot,nky_tot) :: buffer_pt_yt_zc ! full p, full y, constant z
- COMPLEX(dp), DIMENSION(np_loc) :: buffer_pl_cy_zc !local p, constant y, constant z
- COMPLEX(dp), DIMENSION(np_tot) :: buffer_pt_cy_zc ! full p, constant y, constant z
+ INTEGER, INTENT(IN) :: na_tot,np_loc,np_tot,nj_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot
+ COMPLEX(dp), DIMENSION(na_tot,np_loc,nj_tot,nky_loc,nkx_tot,nz_loc), INTENT(IN) :: field_loc
+ COMPLEX(dp), DIMENSION(na_tot,np_tot,nj_tot,nky_tot,nkx_tot,nz_tot), INTENT(OUT) :: field_tot
+ COMPLEX(dp), DIMENSION(na_tot,np_tot,nky_tot,nz_loc) :: buffer_pt_yt_zl ! full p, full y, local z
+ COMPLEX(dp), DIMENSION(na_tot,np_tot,nky_tot,nz_tot) :: buffer_pt_yt_zt ! full p, full y, full z
+ COMPLEX(dp), DIMENSION(na_tot,np_tot,nky_loc) :: buffer_pt_yl_zc ! full p, local y, constant z
+ COMPLEX(dp), DIMENSION(na_tot,np_tot,nky_tot) :: buffer_pt_yt_zc ! full p, full y, constant z
+ COMPLEX(dp), DIMENSION(na_tot,np_loc) :: buffer_pl_cy_zc !local p, constant y, constant z
+ COMPLEX(dp), DIMENSION(na_tot,np_tot) :: buffer_pt_cy_zc ! full p, constant y, constant z
INTEGER :: snd_p, snd_y, snd_z, root_p, root_z, root_ky, iy, ix, iz, ij
- snd_p = np_loc ! Number of points to send along p (per z,y)
- snd_y = np_tot*nky_loc ! Number of points to send along y (per z, full p)
- snd_z = np_tot*nky_tot*nz_loc ! Number of points to send along z (full y, full p)
+ snd_p = na_tot*np_loc ! Number of points to send along p (per z,y)
+ snd_y = na_tot*np_tot*nky_loc ! Number of points to send along y (per z, full p)
+ snd_z = na_tot*np_tot*nky_tot*nz_loc ! Number of points to send along z (full y, full p)
root_p = 0; root_z = 0; root_ky = 0
j: DO ij = 1,nj_tot
x: DO ix = 1,nkx_tot
z: DO iz = 1,nz_loc
y: DO iy = 1,nky_loc
! fill a buffer to contain a slice of p data at constant j, ky, kx and z
- buffer_pl_cy_zc(1:np_loc) = field_loc(1:np_loc,ij,iy,ix,iz)
+ buffer_pl_cy_zc(1:na_tot,1:np_loc) = field_loc(1:na_tot,1:np_loc,ij,iy,ix,iz)
CALL MPI_GATHERV(buffer_pl_cy_zc, snd_p, MPI_DOUBLE_COMPLEX, &
buffer_pt_cy_zc, rcv_p, dsp_p, MPI_DOUBLE_COMPLEX, &
root_p, comm_p, ierr)
- buffer_pt_yl_zc(1:np_tot,iy) = buffer_pt_cy_zc(1:np_tot)
+ buffer_pt_yl_zc(1:na_tot,1:np_tot,iy) = buffer_pt_cy_zc(1:na_tot,1:np_tot)
ENDDO y
! send the full line on p contained by root_p
IF(rank_p .EQ. 0) THEN
CALL MPI_GATHERV(buffer_pt_yl_zc, snd_y, MPI_DOUBLE_COMPLEX, &
buffer_pt_yt_zc, rcv_yp, dsp_yp, MPI_DOUBLE_COMPLEX, &
root_ky, comm_ky, ierr)
- buffer_pt_yt_zl(1:np_tot,1:nky_tot,iz) = buffer_pt_yt_zc(1:np_tot,1:nky_tot)
+ buffer_pt_yt_zl(1:na_tot,1:np_tot,1:nky_tot,iz) = buffer_pt_yt_zc(1:na_tot,1:np_tot,1:nky_tot)
ENDIF
ENDDO z
! send the full line on y contained by root_kyas
@@ -319,7 +319,7 @@ CONTAINS
ENDIF
! ID 0 (the one who ouptut) rebuild the whole array
IF(my_id .EQ. 0) &
- field_tot(1:np_tot,ij,1:nky_tot,ix,1:nz_tot) = buffer_pt_yt_zt(1:np_tot,1:nky_tot,1:nz_tot)
+ field_tot(1:na_tot,1:np_tot,ij,1:nky_tot,ix,1:nz_tot) = buffer_pt_yt_zt(1:na_tot,1:np_tot,1:nky_tot,1:nz_tot)
ENDDO x
ENDDO j
END SUBROUTINE gather_pjxyz
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 001993a6bda953ab2347e0178e4a4e4bcc04facd..0f88cc6b42762ca560a862295605c45539dfb13a 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -1,509 +1,516 @@
MODULE processing
- USE prec_const, ONLY: dp, imagu, SQRT2, SQRT3
- USE grid, ONLY: &
- local_na, local_np, local_nj, local_nky, local_nkx, local_nz, Ngz,Ngj,Ngp, &
- parray,pmax,ip0, iodd, ieven,&
- CONTAINSp0,ip1,CONTAINSp1,ip2,CONTAINSp2,ip3,CONTAINSp3,&
- jarray,jmax,ij0, dmax,&
- kyarray, AA_y,&
- kxarray, AA_x,&
- zarray, deltaz, ieven, iodd, inv_deltaz
- USE fields, ONLY: moments, phi, psi
- USE array, ONLY : kernel, nadiab_moments, &
- ddz_napj, ddzND_Napj, interp_napj,&
- dens, upar, uper, Tpar, Tper, temp
- USE geometry, ONLY: Jacobian, iInt_Jacobian
- USE time_integration, ONLY: updatetlevel
- USE calculus, ONLY: simpson_rule_z, grad_z, grad_z2, grad_z4, interp_z
- USE model, ONLY: EM, CLOS, beta, HDz_h
- USE species, ONLY: tau,q_tau,q_o_sqrt_tau_sigma,sqrt_tau_o_sigma
- USE basic, ONLY: t0_process, t1_process, tc_process
- USE parallel, ONLY: num_procs_ky, rank_ky, comm_ky
- USE mpi
- implicit none
+ USE prec_const, ONLY: dp, imagu, SQRT2, SQRT3
+ USE grid, ONLY: &
+ local_na, local_np, local_nj, local_nky, local_nkx, local_nz, Ngz,Ngj,Ngp, &
+ parray,pmax,ip0, iodd, ieven,&
+ CONTAINSp0,ip1,CONTAINSp1,ip2,CONTAINSp2,ip3,CONTAINSp3,&
+ jarray,jmax,ij0, dmax,&
+ kyarray, AA_y,&
+ kxarray, AA_x,&
+ zarray, deltaz, ieven, iodd, inv_deltaz
+ USE fields, ONLY: moments, phi, psi
+ USE array, ONLY : kernel, nadiab_moments, &
+ ddz_napj, ddzND_Napj, interp_napj,&
+ dens, upar, uper, Tpar, Tper, temp
+ USE geometry, ONLY: Jacobian, iInt_Jacobian
+ USE time_integration, ONLY: updatetlevel
+ USE calculus, ONLY: simpson_rule_z, grad_z, grad_z2, grad_z4, interp_z
+ USE model, ONLY: EM, CLOS, beta, HDz_h
+ USE species, ONLY: tau,q_tau,q_o_sqrt_tau_sigma,sqrt_tau_o_sigma
+ USE basic, ONLY: t0_process, t1_process, tc_process
+ USE parallel, ONLY: num_procs_ky, rank_ky, comm_ky
+ USE mpi
+ implicit none
- REAL(dp), PUBLIC, ALLOCATABLE, DIMENSION(:), PROTECTED :: pflux_x, gflux_x
- REAL(dp), PUBLIC, ALLOCATABLE, DIMENSION(:), PROTECTED :: hflux_x
- INTEGER :: ierr
- PUBLIC :: init_process
- PUBLIC :: compute_nadiab_moments_z_gradients_and_interp
- PUBLIC :: compute_density, compute_upar, compute_uperp
- PUBLIC :: compute_Tpar, compute_Tperp, compute_fluid_moments
- PUBLIC :: compute_radial_transport
- PUBLIC :: compute_radial_heatflux
- PUBLIC :: compute_Napjz_spectrum
+ REAL(dp), PUBLIC, ALLOCATABLE, DIMENSION(:), PROTECTED :: pflux_x, gflux_x
+ REAL(dp), PUBLIC, ALLOCATABLE, DIMENSION(:), PROTECTED :: hflux_x
+ INTEGER :: ierr
+ PUBLIC :: init_process
+ PUBLIC :: compute_nadiab_moments_z_gradients_and_interp
+ PUBLIC :: compute_density, compute_upar, compute_uperp
+ PUBLIC :: compute_Tpar, compute_Tperp, compute_fluid_moments
+ PUBLIC :: compute_radial_transport
+ PUBLIC :: compute_radial_heatflux
+ PUBLIC :: compute_Napjz_spectrum
CONTAINS
-SUBROUTINE init_process
- USE grid, ONLY: local_na
- IMPLICIT NONE
- ALLOCATE( pflux_x(local_na))
- ALLOCATE( gflux_x(local_na))
- ALLOCATE( hflux_x(local_na))
-END SUBROUTINE init_process
+ SUBROUTINE init_process
+ USE grid, ONLY: local_na
+ IMPLICIT NONE
+ ALLOCATE( pflux_x(local_na))
+ ALLOCATE( gflux_x(local_na))
+ ALLOCATE( hflux_x(local_na))
+ END SUBROUTINE init_process
! 1D diagnostic to compute the average radial particle transport <n_a v_ExB_x>_xyz
-SUBROUTINE compute_radial_transport
- IMPLICIT NONE
- COMPLEX(dp) :: pflux_local, gflux_local, integral
- REAL(dp) :: buffer(2)
- INTEGER :: i_, root, iky, ikx, ia, iz, in
- COMPLEX(dp), DIMENSION(local_nz+Ngz) :: integrant
- DO ia = 1,local_na
- pflux_local = 0._dp ! particle flux
- gflux_local = 0._dp ! gyrocenter flux
- integrant = 0._dp ! auxiliary variable for z integration
- !!---------- Gyro center flux (drift kinetic) ------------
- ! Electrostatic part
- IF(CONTAINSp0) THEN
- DO iz = 1,local_nz+ngz ! we include ghost for integration
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- integrant(iz) = integrant(iz) &
- +Jacobian(iz,ieven)*moments(ia,ip0,ij0,iky,ikx,iz,updatetlevel)&
- *imagu*kyarray(iky)*CONJG(phi(iky,ikx,iz))
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- ! Electromagnetic part
- IF( EM .AND. CONTAINSp1 ) THEN
- DO iz = 1,local_nz+ngz
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- integrant(iz) = integrant(iz)&
- -Jacobian(iz,iodd)*sqrt_tau_o_sigma(ia)*moments(ia,ip1,ij0,iky,ikx,iz,updatetlevel)&
- *imagu*kyarray(iky)*CONJG(psi(iky,ikx,iz))
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- ! Integrate over z
- call simpson_rule_z(local_nz,deltaz,integrant,integral)
- ! Get process local gyrocenter flux with a factor two to account for the negative ky modes
- gflux_local = 2._dp*integral*iInt_Jacobian
-
- !
- integrant = 0._dp ! reset auxiliary variable
- !!---------- Particle flux (gyrokinetic) ------------
- ! Electrostatic part
- IF(CONTAINSp0) THEN
- DO iz = 1,local_nz+ngz
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO in = 1+ngj/2, local_nj+ngj/2 ! only interior points
- integrant(iz) = integrant(iz)+ &
- Jacobian(iz,ieven)*moments(ia,ip0,in,iky,ikx,iz,updatetlevel)&
- *imagu*kyarray(iky)*kernel(ia,in,iky,ikx,iz,ieven)*CONJG(phi(iky,ikx,iz))
- ENDDO
+ SUBROUTINE compute_radial_transport
+ IMPLICIT NONE
+ COMPLEX(dp) :: pflux_local, gflux_local, integral
+ REAL(dp) :: buffer(2)
+ INTEGER :: i_, root, iky, ikx, ia, iz, in, izi
+ COMPLEX(dp), DIMENSION(local_nz) :: integrant
+ DO ia = 1,local_na
+ pflux_local = 0._dp ! particle flux
+ gflux_local = 0._dp ! gyrocenter flux
+ integrant = 0._dp ! auxiliary variable for z integration
+ !!---------- Gyro center flux (drift kinetic) ------------
+ ! Electrostatic part
+ IF(CONTAINSp0) THEN
+ DO iz = 1,local_nz
+ izi = iz + ngz/2 !interior index for ghosted arrays
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ integrant(iz) = integrant(iz) &
+ +Jacobian(izi,ieven)*moments(ia,ip0,ij0,iky,ikx,izi,updatetlevel)&
+ *imagu*kyarray(iky)*CONJG(phi(iky,ikx,izi))
+ ENDDO
+ ENDDO
ENDDO
- ENDDO
- ENDDO
- ENDIF
- ! Electromagnetic part
- IF( EM .AND. CONTAINSp1 ) THEN
- DO iz = 1,local_nz+ngz
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO in = 1+ngj/2, local_nj+ngj/2 ! only interior points
- integrant(iz) = integrant(iz) - &
- Jacobian(iz,iodd)*sqrt_tau_o_sigma(ia)*moments(ia,ip1,in,iky,ikx,iz,updatetlevel)&
- *imagu*kyarray(iky)*kernel(ia,in,iky,ikx,iz,iodd)*CONJG(psi(iky,ikx,iz))
+ ENDIF
+ ! Electromagnetic part
+ IF( EM .AND. CONTAINSp1 ) THEN
+ DO iz = 1,local_nz ! we take interior points only
+ izi = iz + ngz/2 !interior index for ghosted arrays
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ integrant(iz) = integrant(iz)&
+ -Jacobian(izi,iodd)*sqrt_tau_o_sigma(ia)*moments(ia,ip1,ij0,iky,ikx,izi,updatetlevel)&
+ *imagu*kyarray(iky)*CONJG(psi(iky,ikx,izi))
+ ENDDO
+ ENDDO
ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- ! Integrate over z
- call simpson_rule_z(local_nz,deltaz,integrant,integral)
- ! Get process local particle flux with a factor two to account for the negative ky modes
- pflux_local = 2._dp*integral*iInt_Jacobian
+ ENDIF
+ ! Integrate over z
+ call simpson_rule_z(local_nz,deltaz,integrant,integral)
+ ! Get process local gyrocenter flux with a factor two to account for the negative ky modes
+ gflux_local = 2._dp*integral*iInt_Jacobian
- !!!!---------- Sum over all processes ----------
- buffer(1) = REAL(gflux_local)
- buffer(2) = REAL(pflux_local)
- root = 0
- !Gather manually among the rank_p=0 processes and perform the sum
- gflux_x(ia) = 0
- pflux_x(ia) = 0
- 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, 2 , MPI_DOUBLE_PRECISION, root, 1234, 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, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
- gflux_x(ia) = gflux_x(ia) + buffer(1)
- pflux_x(ia) = pflux_x(ia) + buffer(2)
+ !
+ integrant = 0._dp ! reset auxiliary variable
+ !!---------- Particle flux (gyrokinetic) ------------
+ ! Electrostatic part
+ IF(CONTAINSp0) THEN
+ DO iz = 1,local_nz ! we take interior points only
+ izi = iz + ngz/2 !interior index for ghosted arrays
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO in = 1+ngj/2, local_nj+ngj/2 ! only interior points
+ integrant(iz) = integrant(iz)+ &
+ Jacobian(izi,ieven)*moments(ia,ip0,in,iky,ikx,izi,updatetlevel)&
+ *imagu*kyarray(iky)*kernel(ia,in,iky,ikx,izi,ieven)*CONJG(phi(iky,ikx,izi))
+ ENDDO
+ ENDDO
+ ENDDO
ENDDO
- ENDIF
- ELSE
- gflux_x(ia) = gflux_local
- pflux_x(ia) = pflux_local
- ENDIF
- ENDDO
- ! if(my_id .eq. 0) write(*,*) 'pflux_ri = ',pflux_ri
-END SUBROUTINE compute_radial_transport
+ ENDIF
+ ! Electromagnetic part
+ IF( EM .AND. CONTAINSp1 ) THEN
+ DO iz = 1,local_nz ! we take interior points only
+ izi = iz + ngz/2 !interior index for ghosted arrays
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO in = 1+ngj/2, local_nj+ngj/2 ! only interior points
+ integrant(iz) = integrant(iz) - &
+ Jacobian(izi,iodd)*sqrt_tau_o_sigma(ia)*moments(ia,ip1,in,iky,ikx,izi,updatetlevel)&
+ *imagu*kyarray(iky)*kernel(ia,in,iky,ikx,izi,iodd)*CONJG(psi(iky,ikx,izi))
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ! Integrate over z
+ call simpson_rule_z(local_nz,deltaz,integrant,integral)
+ ! Get process local particle flux with a factor two to account for the negative ky modes
+ pflux_local = 2._dp*integral*iInt_Jacobian
-! 1D diagnostic to compute the average radial particle transport <T_i v_ExB_x>_xyz
-SUBROUTINE compute_radial_heatflux
- IMPLICIT NONE
- COMPLEX(dp) :: hflux_local, integral
- REAL(dp) :: buffer(2), n_dp
- INTEGER :: i_, root, in, ia, iky, ikx, iz
- COMPLEX(dp), DIMENSION(local_nz+ngz) :: integrant ! charge density q_a n_a
- DO ia = 1,local_na
- hflux_local = 0._dp ! heat flux
- integrant = 0._dp ! z integration auxiliary variable
- !!----------------ELECTROSTATIC CONTRIBUTION---------------------------
- IF(CONTAINSp0 .AND. CONTAINSp2) THEN
- ! Loop to compute gamma_kx = sum_ky sum_j -i k_z Kernel_j Na00 * phi
- DO iz = 1,local_nz+ngz
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO in = 1+ngj/2, local_nj+ngj/2 ! only interior points
- n_dp = jarray(in)
- integrant(iz) = integrant(iz) &
- +Jacobian(iz,ieven)*tau(ia)*imagu*kyarray(iky)*CONJG(phi(iky,ikx,iz))&
- *kernel(ia,in,iky,ikx,iz,ieven)*(&
- 0.5_dp*SQRT2*moments(ia,ip2,in ,iky,ikx,iz,updatetlevel)&
- +(2._dp*n_dp + 1.5_dp)*moments(ia,ip0,in ,iky,ikx,iz,updatetlevel)&
- -(n_dp+1._dp)*moments(ia,ip0,in+1,iky,ikx,iz,updatetlevel)&
- -n_dp*moments(ia,ip0,in-1,iky,ikx,iz,updatetlevel))
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- IF(EM .AND. CONTAINSp1 .AND. CONTAINSp3) THEN
- !!----------------ELECTROMAGNETIC CONTRIBUTION--------------------
- DO iz = 1,local_nz+ngz
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO in = 1+ngj/2, local_nj+ngj/2 ! only interior points
- n_dp = jarray(in)
- integrant(iz) = integrant(iz) &
- +Jacobian(iz,iodd)*tau(ia)*sqrt_tau_o_sigma(ia)*imagu*kyarray(iky)*CONJG(psi(iky,ikx,iz))&
- *kernel(ia,in,iky,ikx,iz,iodd)*(&
- 0.5_dp*SQRT2*SQRT3*moments(ia,ip3,in ,iky,ikx,iz,updatetlevel)&
- +1.5_dp*moments(ia,ip1,in ,iky,ikx,iz,updatetlevel)&
- +(2._dp*n_dp+1._dp)*moments(ia,ip1,in ,iky,ikx,iz,updatetlevel)&
- -(n_dp+1._dp)*moments(ia,ip1,in+1,iky,ikx,iz,updatetlevel)&
- -n_dp*moments(ia,ip1,in-1,iky,ikx,iz,updatetlevel))
+ !!!!---------- Sum over all processes ----------
+ buffer(1) = REAL(gflux_local,dp)
+ buffer(2) = REAL(pflux_local,dp)
+ root = 0
+ !Gather manually among the rank_p=0 processes and perform the sum
+ gflux_x(ia) = 0
+ pflux_x(ia) = 0
+ 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, 2 , MPI_DOUBLE_PRECISION, root, 1234, 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, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
+ gflux_x(ia) = gflux_x(ia) + buffer(1)
+ pflux_x(ia) = pflux_x(ia) + buffer(2)
+ ENDDO
+ ENDIF
+ ELSE
+ gflux_x(ia) = gflux_local
+ pflux_x(ia) = pflux_local
+ ENDIF
ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- ! Add polarisation contribution
- ! integrant(iz) = integrant(iz) + tau_i*imagu*ky_&
- ! *CONJG(phi(iky,ikx,iz))*phi(iky,ikx,iz) * HF_phi_correction_operator(iky,ikx,iz)
- ! Integrate over z
- call simpson_rule_z(local_nz,deltaz,integrant,integral)
- ! Double it for taking into account the other half plane
- hflux_local = 2._dp*integral*iInt_Jacobian
- buffer(2) = REAL(hflux_local)
- root = 0
- !Gather manually among the rank_p=0 processes and perform the sum
- hflux_x(ia) = 0
- 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, 2 , MPI_DOUBLE_PRECISION, root, 1234, 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, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
+ ! if(my_id .eq. 0) write(*,*) 'pflux_ri = ',pflux_ri
+ END SUBROUTINE compute_radial_transport
+
+! 1D diagnostic to compute the average radial particle transport <T_i v_ExB_x>_xyz
+ SUBROUTINE compute_radial_heatflux
+ IMPLICIT NONE
+ COMPLEX(dp) :: hflux_local, integral
+ REAL(dp) :: buffer(2), n_dp
+ INTEGER :: i_, root, in, ia, iky, ikx, iz, izi
+ COMPLEX(dp), DIMENSION(local_nz) :: integrant ! charge density q_a n_a
+ DO ia = 1,local_na
+ hflux_local = 0._dp ! heat flux
+ integrant = 0._dp ! z integration auxiliary variable
+ !!----------------ELECTROSTATIC CONTRIBUTION---------------------------
+ IF(CONTAINSp0 .AND. CONTAINSp2) THEN
+ ! Loop to compute gamma_kx = sum_ky sum_j -i k_z Kernel_j Na00 * phi
+ DO iz = 1,local_nz ! we take interior points only
+ izi = iz + ngz/2 !interior index for ghosted arrays
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO in = 1+ngj/2, local_nj+ngj/2 ! only interior points
+ n_dp = jarray(in)
+ integrant(iz) = integrant(iz) &
+ +Jacobian(izi,ieven)*tau(ia)*imagu*kyarray(iky)*CONJG(phi(iky,ikx,izi))&
+ *kernel(ia,in,iky,ikx,izi,ieven)*(&
+ 0.5_dp*SQRT2*moments(ia,ip2,in ,iky,ikx,izi,updatetlevel)&
+ +(2._dp*n_dp + 1.5_dp)*moments(ia,ip0,in ,iky,ikx,izi,updatetlevel)&
+ -(n_dp+1._dp)*moments(ia,ip0,in+1,iky,ikx,izi,updatetlevel)&
+ -n_dp*moments(ia,ip0,in-1,iky,ikx,izi,updatetlevel))
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ IF(EM .AND. CONTAINSp1 .AND. CONTAINSp3) THEN
+ !!----------------ELECTROMAGNETIC CONTRIBUTION--------------------
+ DO iz = 1,local_nz
+ izi = iz + ngz/2 !interior index for ghosted arrays
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO in = 1+ngj/2, local_nj+ngj/2 ! only interior points
+ n_dp = jarray(in)
+ integrant(iz) = integrant(iz) &
+ +Jacobian(izi,iodd)*tau(ia)*sqrt_tau_o_sigma(ia)*imagu*kyarray(iky)*CONJG(psi(iky,ikx,izi))&
+ *kernel(ia,in,iky,ikx,izi,iodd)*(&
+ 0.5_dp*SQRT2*SQRT3*moments(ia,ip3,in ,iky,ikx,izi,updatetlevel)&
+ +1.5_dp*moments(ia,ip1,in ,iky,ikx,izi,updatetlevel)&
+ +(2._dp*n_dp+1._dp)*moments(ia,ip1,in ,iky,ikx,izi,updatetlevel)&
+ -(n_dp+1._dp)*moments(ia,ip1,in+1,iky,ikx,izi,updatetlevel)&
+ -n_dp*moments(ia,ip1,in-1,iky,ikx,izi,updatetlevel))
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ! Add polarisation contribution
+ ! integrant(iz) = integrant(iz) + tau_i*imagu*ky_&
+ ! *CONJG(phi(iky,ikx,iz))*phi(iky,ikx,iz) * HF_phi_correction_operator(iky,ikx,iz)
+ ! Integrate over z
+ call simpson_rule_z(local_nz,deltaz,integrant,integral)
+ ! Double it for taking into account the other half plane
+ hflux_local = 2._dp*integral*iInt_Jacobian
+ buffer(2) = REAL(hflux_local,dp)
+ root = 0
+ !Gather manually among the rank_p=0 processes and perform the sum
+ hflux_x(ia) = 0
+ 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, 2 , MPI_DOUBLE_PRECISION, root, 1234, 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, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
hflux_x(ia) = hflux_x(ia) + buffer(2)
- ENDDO
- ENDIF
- ELSE
- hflux_x(ia) = hflux_local
- ENDIF
- ENDDO
-END SUBROUTINE compute_radial_heatflux
+ ENDDO
+ ENDIF
+ ELSE
+ hflux_x(ia) = hflux_local
+ ENDIF
+ ENDDO
+ END SUBROUTINE compute_radial_heatflux
-SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
- ! evaluate the non-adiabatique ion moments
- !
- ! n_{pi} = N^{pj} + kernel(j) /tau_i phi delta_p0
- !
- IMPLICIT NONE
- INTEGER :: eo, p_int, j_int, ia,ip,ij,iky,ikx,iz
- COMPLEX(dp), DIMENSION(local_nz+ngz) :: f_in
- COMPLEX(dp), DIMENSION(local_nz) :: f_out
- CALL cpu_time(t0_process)
+ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
+ ! evaluate the non-adiabatique ion moments
+ !
+ ! n_{pi} = N^{pj} + kernel(j) /tau_i phi delta_p0
+ !
+ IMPLICIT NONE
+ INTEGER :: eo, p_int, j_int, ia,ip,ij,iky,ikx,iz
+ COMPLEX(dp), DIMENSION(local_nz+ngz) :: f_in
+ COMPLEX(dp), DIMENSION(local_nz) :: f_out
+ CALL cpu_time(t0_process)
- !non adiab moments
- DO iz=1,local_nz+ngz
- DO ikx=1,local_nkx
- DO iky=1,local_nky
- DO ij=1,local_nj+ngj
- DO ip=1,local_np+ngp
- DO ia = 1,local_na
- IF(parray(ip) .EQ. 0) THEN
- nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel) &
- + q_tau(ia)*kernel(ia,ij,iky,ikx,iz,ieven)*phi(iky,ikx,iz)
- ELSEIF( (parray(ip) .EQ. 1) .AND. (beta .GT. 0) ) THEN
- nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel) &
- - q_o_sqrt_tau_sigma(ia)*kernel(ia,ij,iky,ikx,iz,ieven)*psi(iky,ikx,iz)
- ELSE
- nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel)
- ENDIF
+ !non adiab moments
+ DO iz=1,local_nz+ngz
+ DO ikx=1,local_nkx
+ DO iky=1,local_nky
+ DO ij=1,local_nj+ngj
+ DO ip=1,local_np+ngp
+ DO ia = 1,local_na
+ IF(parray(ip) .EQ. 0) THEN
+ nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel) &
+ + q_tau(ia)*kernel(ia,ij,iky,ikx,iz,ieven)*phi(iky,ikx,iz)
+ ELSEIF( (parray(ip) .EQ. 1) .AND. (beta .GT. 0) ) THEN
+ nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel) &
+ - q_o_sqrt_tau_sigma(ia)*kernel(ia,ij,iky,ikx,iz,ieven)*psi(iky,ikx,iz)
+ ELSE
+ nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel)
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
ENDDO
- ENDDO
- ENDDO
+ ENDDO
ENDDO
- ENDDO
- ENDDO
- !! Ensure to kill the moments too high if closue option is set to 1
- IF(CLOS .EQ. 1) THEN
- DO ij=1,local_nj+ngj
- j_int = jarray(ij)
- DO ip=1,local_np+ngp
- p_int = parray(ip)
- DO ia = 1,local_na
- IF(p_int+2*j_int .GT. dmax) &
- nadiab_moments(ia,ip,ij,:,:,:) = 0._dp
- ENDDO
- ENDDO
- ENDDO
- ENDIF
+ !! Ensure to kill the moments too high if closue option is set to 1
+ IF(CLOS .EQ. 1) THEN
+ DO ij=1,local_nj+ngj
+ j_int = jarray(ij)
+ DO ip=1,local_np+ngp
+ p_int = parray(ip)
+ DO ia = 1,local_na
+ IF(p_int+2*j_int .GT. dmax) &
+ nadiab_moments(ia,ip,ij,:,:,:) = 0._dp
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
- !------------- INTERP AND GRADIENTS ALONG Z ----------------------------------
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO ij = 1,local_nj+ngj
- DO ip = 1,local_np+ngp
- DO ia = 1,local_na
- p_int = parray(ip)
- eo = MODULO(p_int,2)+1 ! Indicates if we are on even or odd z grid
- ! Compute z first derivative
- f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
- CALL grad_z(eo,local_nz,ngz,inv_deltaz,f_in,f_out)
- ddz_napj(ia,ip,ij,iky,ikx,:) = f_out
- ! Parallel numerical diffusion
- IF (HDz_h) THEN
- f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
- CALL grad_z4(local_nz,ngz,inv_deltaz,f_in,f_out)
- ddzND_Napj(ia,ip,ij,iky,ikx,:) = f_out
- ELSE
- f_in = moments(ia,ip,ij,iky,ikx,:,updatetlevel)
- CALL grad_z4(local_nz,ngz,inv_deltaz,f_in,f_out)
- ddzND_Napj(ia,ip,ij,iky,ikx,:) = f_out
- ENDIF
- ! Compute even odd grids interpolation
- f_in = nadiab_moments(ia,ip,ij,iky,ikx,1:local_nz+ngz)
- CALL interp_z(eo,local_nz,ngz,f_in,f_out)
- interp_napj(ia,ip,ij,iky,ikx,1:local_nz) = f_out
+ !------------- INTERP AND GRADIENTS ALONG Z ----------------------------------
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO ij = 1,local_nj+ngj
+ DO ip = 1,local_np+ngp
+ DO ia = 1,local_na
+ p_int = parray(ip)
+ eo = MODULO(p_int,2)+1 ! Indicates if we are on even or odd z grid
+ ! Compute z first derivative
+ f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
+ CALL grad_z(eo,local_nz,ngz,inv_deltaz,f_in,f_out)
+ ddz_napj(ia,ip,ij,iky,ikx,:) = f_out
+ ! Parallel numerical diffusion
+ IF (HDz_h) THEN
+ f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
+ CALL grad_z4(local_nz,ngz,inv_deltaz,f_in,f_out)
+ ddzND_Napj(ia,ip,ij,iky,ikx,:) = f_out
+ ELSE
+ f_in = moments(ia,ip,ij,iky,ikx,:,updatetlevel)
+ CALL grad_z4(local_nz,ngz,inv_deltaz,f_in,f_out)
+ ddzND_Napj(ia,ip,ij,iky,ikx,:) = f_out
+ ENDIF
+ ! Compute even odd grids interpolation
+ f_in = nadiab_moments(ia,ip,ij,iky,ikx,1:local_nz+ngz)
+ CALL interp_z(eo,local_nz,ngz,f_in,f_out)
+ interp_napj(ia,ip,ij,iky,ikx,1:local_nz) = f_out
+ ENDDO
+ ENDDO
ENDDO
- ENDDO
- ENDDO
+ ENDDO
ENDDO
- ENDDO
- ! Phi parallel gradient (not implemented fully, should be negligible)
- ! DO ikx = 1,local_nkx
- ! DO iky = 1,local_nky
- ! CALL grad_z(0,phi(iky,ikx,1:local_nz+ngz), ddz_phi(iky,ikx,1:local_nz))
- ! ENDDO
- ! ENDDO
+ ! Phi parallel gradient (not implemented fully, should be negligible)
+ ! DO ikx = 1,local_nkx
+ ! DO iky = 1,local_nky
+ ! CALL grad_z(0,phi(iky,ikx,1:local_nz+ngz), ddz_phi(iky,ikx,1:local_nz))
+ ! ENDDO
+ ! ENDDO
- ! Execution time end
- CALL cpu_time(t1_process)
- tc_process = tc_process + (t1_process - t0_process)
-END SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
+ ! Execution time end
+ CALL cpu_time(t1_process)
+ tc_process = tc_process + (t1_process - t0_process)
+ END SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
-SUBROUTINE compute_Napjz_spectrum
- USE fields, ONLY : moments
- USE array, ONLY : Napjz
- USE time_integration, ONLY : updatetlevel
- IMPLICIT NONE
- REAL(dp), DIMENSION(local_np,local_nj,local_nz) :: local_sum,global_sum, buffer
- INTEGER :: i_, root, count, ia, ip, ij, iky, ikx, iz
- root = 0
- DO ia=1,local_na
- ! z-moment spectrum
- ! build local sum
- local_sum = 0._dp
- DO iz = 1,local_nz
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO ij = 1,local_nj
- DO ip = 1,local_np
- local_sum(ip,ij,iz) = local_sum(ip,ij,iz) + &
- (moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel) &
- * CONJG(moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel)))
+ SUBROUTINE compute_Napjz_spectrum
+ USE fields, ONLY : moments
+ USE array, ONLY : Napjz
+ USE time_integration, ONLY : updatetlevel
+ IMPLICIT NONE
+ REAL(dp), DIMENSION(local_np,local_nj,local_nz) :: local_sum,global_sum, buffer
+ INTEGER :: i_, root, count, ia, ip, ij, iky, ikx, iz
+ root = 0
+ DO ia=1,local_na
+ ! z-moment spectrum
+ ! build local sum
+ local_sum = 0._dp
+ DO iz = 1,local_nz
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO ij = 1,local_nj
+ DO ip = 1,local_np
+ local_sum(ip,ij,iz) = local_sum(ip,ij,iz) + &
+ (moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel) &
+ * CONJG(moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel)))
+ ENDDO
+ ENDDO
+ ENDDO
ENDDO
- ENDDO
- ENDDO
+ ENDDO
+ ! sum reduction
+ buffer = local_sum
+ global_sum = 0._dp
+ count = local_np*local_nj*local_nz
+ 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, 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, count , MPI_DOUBLE_PRECISION, i_, 5678, comm_ky, MPI_STATUS_IGNORE, ierr)
+ global_sum = global_sum + buffer
+ ENDDO
+ ENDIF
+ ELSE
+ global_sum = local_sum
+ ENDIF
+ Napjz(ia,:,:,:) = global_sum
ENDDO
- ENDDO
- ! sum reduction
- buffer = local_sum
- global_sum = 0._dp
- count = local_np*local_nj*local_nz
- 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, 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, count , MPI_DOUBLE_PRECISION, i_, 5678, comm_ky, MPI_STATUS_IGNORE, ierr)
- global_sum = global_sum + buffer
- ENDDO
- ENDIF
- ELSE
- global_sum = local_sum
- ENDIF
- Napjz(ia,:,:,:) = global_sum
- ENDDO
-END SUBROUTINE compute_Napjz_spectrum
+ END SUBROUTINE compute_Napjz_spectrum
!_____________________________________________________________________________!
!!!!! FLUID MOMENTS COMPUTATIONS !!!!!
! Compute the 2D particle density for electron and ions (sum over Laguerre)
-SUBROUTINE compute_density
- IMPLICIT NONE
- COMPLEX(dp) :: dens_
- INTEGER :: ia, iz, iky, ikx, ij
- DO ia=1,local_na
- IF ( CONTAINSp0 ) THEN
- ! Loop to compute dens_i = sum_j kernel_j Ni0j at each k
- DO iz = 1,local_nz
- DO iky = 1,local_nky
- DO ikx = 1,local_nkx
- dens_ = 0._dp
- DO ij = 1, local_nj
- dens_ = dens_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven) * moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
+ SUBROUTINE compute_density
+ IMPLICIT NONE
+ COMPLEX(dp) :: dens_
+ INTEGER :: ia, iz, iky, ikx, ij
+ DO ia=1,local_na
+ IF ( CONTAINSp0 ) THEN
+ ! Loop to compute dens_i = sum_j kernel_j Ni0j at each k
+ DO iz = 1,local_nz
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ dens_ = 0._dp
+ DO ij = 1, local_nj
+ dens_ = dens_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven) * moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
+ ENDDO
+ dens(ia,iky,ikx,iz) = dens_
+ ENDDO
+ ENDDO
ENDDO
- dens(ia,iky,ikx,iz) = dens_
- ENDDO
- ENDDO
+ ENDIF
ENDDO
- ENDIF
- ENDDO
-END SUBROUTINE compute_density
+ END SUBROUTINE compute_density
! Compute the 2D particle fluid perp velocity for electron and ions (sum over Laguerre)
-SUBROUTINE compute_uperp
- IMPLICIT NONE
- COMPLEX(dp) :: uperp_
- INTEGER :: ia, iz, iky, ikx, ij
- DO ia=1,local_na
- IF ( CONTAINSp0 ) THEN
- DO iz = 1,local_nz
- DO iky = 1,local_nky
- DO ikx = 1,local_nkx
- uperp_ = 0._dp
- DO ij = 1, local_nj
- uperp_ = uperp_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven) *&
- 0.5_dp*(moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel) - moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
+ SUBROUTINE compute_uperp
+ IMPLICIT NONE
+ COMPLEX(dp) :: uperp_
+ INTEGER :: ia, iz, iky, ikx, ij
+ DO ia=1,local_na
+ IF ( CONTAINSp0 ) THEN
+ DO iz = 1,local_nz
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ uperp_ = 0._dp
+ DO ij = 1, local_nj
+ uperp_ = uperp_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven) *&
+ 0.5_dp*(moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
+ -moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
+ ENDDO
+ uper(ia,iky,ikx,iz) = uperp_
+ ENDDO
ENDDO
- uper(ia,iky,ikx,iz) = uperp_
ENDDO
- ENDDO
- ENDDO
- ENDIF
- ENDDO
-END SUBROUTINE compute_uperp
+ ENDIF
+ ENDDO
+ END SUBROUTINE compute_uperp
! Compute the 2D particle fluid par velocity for electron and ions (sum over Laguerre)
-SUBROUTINE compute_upar
- IMPLICIT NONE
- INTEGER :: ia, iz, iky, ikx, ij
- COMPLEX(dp) :: upar_
- DO ia=1,local_na
- IF ( CONTAINSp1 ) THEN
- DO iz = 1,local_nz
- DO iky = 1,local_nky
- DO ikx = 1,local_nkx
- upar_ = 0._dp
- DO ij = 1, local_nj
- upar_ = upar_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,iodd)*moments(ia,ip1,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
- ENDDO
- upar(ia,iky,ikx,iz) = upar_
- ENDDO
- ENDDO
+ SUBROUTINE compute_upar
+ IMPLICIT NONE
+ INTEGER :: ia, iz, iky, ikx, ij
+ COMPLEX(dp) :: upar_
+ DO ia=1,local_na
+ IF ( CONTAINSp1 ) THEN
+ DO iz = 1,local_nz
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ upar_ = 0._dp
+ DO ij = 1, local_nj
+ upar_ = upar_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,iodd)*moments(ia,ip1,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
+ ENDDO
+ upar(ia,iky,ikx,iz) = upar_
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
ENDDO
- ENDIF
- ENDDO
-END SUBROUTINE compute_upar
+ END SUBROUTINE compute_upar
! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
-SUBROUTINE compute_tperp
- USE time_integration, ONLY : updatetlevel
- IMPLICIT NONE
- REAL(dp) :: j_dp
- COMPLEX(dp) :: Tperp_
- INTEGER :: ia, iz, iky, ikx, ij
- DO ia=1,local_na
- IF ( CONTAINSp0 .AND. CONTAINSp2 ) THEN
- ! Loop to compute T = 1/3*(Tpar + 2Tperp)
- DO iz = 1,local_nz
- DO iky = 1,local_nky
- DO ikx = 1,local_nkx
- Tperp_ = 0._dp
- DO ij = 1, local_nj
- j_dp = REAL(ij-1,dp)
- Tperp_ = Tperp_ + kernel(ia,ij,iky,ikx,iz,ieven)*&
- ((2_dp*j_dp+1)*moments(ia,ip0,ij +ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
- -j_dp *moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
- -j_dp+1 *moments(ia,ip0,ij+1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
- ENDDO
- Tper(ia,iky,ikx,iz) = Tperp_
+ SUBROUTINE compute_tperp
+ USE time_integration, ONLY : updatetlevel
+ IMPLICIT NONE
+ REAL(dp) :: j_dp
+ COMPLEX(dp) :: Tperp_
+ INTEGER :: ia, iz, iky, ikx, ij
+ DO ia=1,local_na
+ IF ( CONTAINSp0 .AND. CONTAINSp2 ) THEN
+ ! Loop to compute T = 1/3*(Tpar + 2Tperp)
+ DO iz = 1,local_nz
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ Tperp_ = 0._dp
+ DO ij = 1, local_nj
+ j_dp = jarray(ij+ngj/2)
+ Tperp_ = Tperp_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven)*&
+ ((2_dp*j_dp+1)*moments(ia,ip0,ij +ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
+ -j_dp*moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
+ -(j_dp+1)*moments(ia,ip0,ij+1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
+ ENDDO
+ Tper(ia,iky,ikx,iz) = Tperp_
+ ENDDO
+ ENDDO
ENDDO
- ENDDO
- ENDDO
- ENDIF
- ENDDO
-END SUBROUTINE compute_Tperp
+ ENDIF
+ ENDDO
+ END SUBROUTINE compute_Tperp
! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
-SUBROUTINE compute_Tpar
- USE time_integration, ONLY : updatetlevel
- IMPLICIT NONE
- REAL(dp) :: j_dp
- COMPLEX(dp) :: Tpar_
- INTEGER :: ia, iz, iky, ikx, ij
+ SUBROUTINE compute_Tpar
+ USE time_integration, ONLY : updatetlevel
+ IMPLICIT NONE
+ REAL(dp) :: j_dp
+ COMPLEX(dp) :: Tpar_
+ INTEGER :: ia, iz, iky, ikx, ij
- DO ia=1,local_na
- IF ( CONTAINSp0 .AND. CONTAINSp0 ) THEN
- ! Loop to compute T = 1/3*(Tpar + 2Tperp)
- DO iz = 1,local_nz
- DO iky = 1,local_nky
- DO ikx = 1,local_nkx
- Tpar_ = 0._dp
- DO ij = 1, local_nj
- j_dp = REAL(ij-1,dp)
- Tpar_ = Tpar_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven)*&
- (SQRT2 * moments(ia,ip2,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel) &
- + moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
+ DO ia=1,local_na
+ IF ( CONTAINSp0 .AND. CONTAINSp0 ) THEN
+ ! Loop to compute T = 1/3*(Tpar + 2Tperp)
+ DO iz = 1,local_nz
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ Tpar_ = 0._dp
+ DO ij = 1, local_nj
+ j_dp = REAL(ij-1,dp)
+ Tpar_ = Tpar_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven)*&
+ (SQRT2 * moments(ia,ip2,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel) &
+ + moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
+ ENDDO
+ Tpar(ia,iky,ikx,iz) = Tpar_
+ ENDDO
+ ENDDO
ENDDO
- Tpar(ia,iky,ikx,iz) = Tpar_
- ENDDO
- ENDDO
+ ENDIF
ENDDO
- ENDIF
- ENDDO
-END SUBROUTINE compute_Tpar
+ END SUBROUTINE compute_Tpar
! Compute the 2D particle fluid moments for electron and ions (sum over Laguerre)
-SUBROUTINE compute_fluid_moments
- IMPLICIT NONE
- CALL compute_density
- CALL compute_upar
- CALL compute_uperp
- CALL compute_Tpar
- CALL compute_Tperp
- ! Temperature
- temp = (Tpar - 2._dp * Tper)/3._dp - dens
-END SUBROUTINE compute_fluid_moments
+ SUBROUTINE compute_fluid_moments
+ IMPLICIT NONE
+ CALL compute_density
+ CALL compute_upar
+ CALL compute_uperp
+ CALL compute_Tpar
+ CALL compute_Tperp
+ ! Temperature
+ temp = (Tpar - 2._dp * Tper)/3._dp - dens
+ END SUBROUTINE compute_fluid_moments
END MODULE processing
diff --git a/src/solve_EM_fields.F90 b/src/solve_EM_fields.F90
index 5e61c5dd97961fe8a88f4c6497bf11e8657f3619..0b78a8dda08005cc91003aef0e375a191953e1bc 100644
--- a/src/solve_EM_fields.F90
+++ b/src/solve_EM_fields.F90
@@ -39,7 +39,7 @@ CONTAINS
x:DO ikx = 1,local_nky
y:DO iky = 1,local_nkx
!!!!!!!!!!!!!!! Compute particle charge density q_a n_a for each evolved species
- DO iz = 1,local_Nz
+ DO iz = 1,local_nz
rho(iz) = 0._dp
DO in=1,total_nj
DO ia = 1,local_na
@@ -59,7 +59,7 @@ CONTAINS
IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (y-average)
! Prepare integrant for z-average
integrant(iz) = Jacobian(iz+ngz/2,ieven)*rho(iz)*inv_pol_ion(iky,ikx,iz)
- call simpson_rule_z(local_Nz,deltaz,integrant,intf_) ! get the flux averaged phi
+ call simpson_rule_z(local_nz,deltaz,integrant,intf_) ! get the flux averaged phi
fsa_phi = intf_ * iInt_Jacobian !Normalize by 1/int(Jxyz)dz
ENDIF
rho(iz) = rho(iz) + fsa_phi
@@ -68,7 +68,7 @@ CONTAINS
! IF (ADIAB_I) THEN
! ENDIF
!!!!!!!!!!!!!!! Inverting the poisson equation
- DO iz = 1,local_Nz
+ DO iz = 1,local_nz
phi(iky,ikx,iz+ngz/2) = inv_poisson_op(iky,ikx,iz)*rho(iz)
ENDDO
ENDDO y
diff --git a/src/stepon.F90 b/src/stepon.F90
index f88353964791cc6c002823e501b80d25d5ae7e1f..b193d207053493197bd2f58c443c655a8c9e7a7d 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -6,6 +6,7 @@ SUBROUTINE stepon
USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
use mpi, ONLY: MPI_COMM_WORLD
USE time_integration, ONLY: ntimelevel
+ USE prec_const, ONLY: dp
IMPLICIT NONE
INTEGER :: num_step, ierr
@@ -140,7 +141,7 @@ SUBROUTINE stepon
moments(ia,ip,ij,iky0,ikx,iz,:) = CONJG(moments(ia,ip,ij,iky0,total_nkx+2-ikx,iz,:))
END DO
! must be real at origin and top right
- moments(ia,ip,ij, iky0,ikx0,iz,:) = REAL(moments(ia,ip,ij, iky0,ikx0,iz,:))
+ moments(ia,ip,ij, iky0,ikx0,iz,:) = REAL(moments(ia,ip,ij, iky0,ikx0,iz,:),dp)
ENDDO a
ENDDO p
ENDDO j
@@ -150,13 +151,13 @@ SUBROUTINE stepon
phi(iky0,ikx,:) = phi(iky0,total_nkx+2-ikx,:)
END DO
! must be real at origin
- phi(iky0,ikx0,:) = REAL(phi(iky0,ikx0,:))
+ phi(iky0,ikx0,:) = REAL(phi(iky0,ikx0,:),dp)
! Psi
DO ikx=2,total_nkx/2 !symmetry at ky = 0
psi(iky0,ikx,:) = psi(iky0,total_nkx+2-ikx,:)
END DO
! must be real at origin
- psi(iky0,ikx0,:) = REAL(psi(iky0,ikx0,:))
+ psi(iky0,ikx0,:) = REAL(psi(iky0,ikx0,:),dp)
ENDIF
END SUBROUTINE enforce_symmetry
diff --git a/src/tesend.F90 b/src/tesend.F90
index ebbdbf3926002dee4a94fa820865b24c33f6a4b2..6d53392c3e34a0744bddbc2fd03cbc1e4fc07052 100644
--- a/src/tesend.F90
+++ b/src/tesend.F90
@@ -24,7 +24,7 @@ SUBROUTINE tesend
!________________________________________________________________________________
! 2. Test on NRUN
- nlend = step .GT. nrun
+ nlend = step .GE. nrun
IF ( nlend ) THEN
CALL speak('NRUN steps done')
RETURN
diff --git a/testcases/smallest_problem/fort.90 b/testcases/smallest_problem/fort.90
index 40604ef74fccb99dfffbb576a349cf1c46c12ea2..fdb25d1d72464ce8c16d1b3b84fa1e7c089cbb6c 100644
--- a/testcases/smallest_problem/fort.90
+++ b/testcases/smallest_problem/fort.90
@@ -91,8 +91,8 @@
&INITIAL_CON
INIT_OPT = 'mom00'
ACT_ON_MODES = 'donothing'
- init_background = 1.0
- init_noiselvl = 0.0
+ init_background = 0.0
+ init_noiselvl = 1.0
iseed = 42
/
&TIME_INTEGRATION_PAR
diff --git a/testcases/smallest_problem/fort_00.90 b/testcases/smallest_problem/fort_00.90
index 9ea911dec07f785af19e624808d6be6e66ff27ee..d37fc3f1b792dd5e9880db12839564461fe4463f 100644
--- a/testcases/smallest_problem/fort_00.90
+++ b/testcases/smallest_problem/fort_00.90
@@ -73,8 +73,8 @@
&INITIAL_CON
INIT_OPT = 'mom00'
ACT_ON_MODES = 'donothing'
- init_background = 1.0
- init_noiselvl = 0.0
+ init_background = 0.0
+ init_noiselvl = 1.0
iseed = 42
/
&TIME_INTEGRATION_PAR