diff --git a/Makefile b/Makefile
index c0c7dee375864f3df27bed217e115b70b4c72752..145be064282e76d61ffdf3b07e9ad01064d2b011 100644
--- a/Makefile
+++ b/Makefile
@@ -104,7 +104,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
$(OBJDIR)/advance_field_mod.o : src/advance_field_mod.F90 \
$(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/initial_par_mod.o \
$(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o \
- $(OBJDIR)/fields_mod.o
+ $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/advance_field_mod.F90 -o $@
$(OBJDIR)/array_mod.o : src/array_mod.F90 \
@@ -189,8 +189,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ghosts_mod.F90 -o $@
$(OBJDIR)/grid_mod.o : src/grid_mod.F90 \
- $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/model_mod.o \
- $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/grid_mod.F90 -o $@
$(OBJDIR)/inital.o : src/inital.F90 \
diff --git a/src/advance_field_mod.F90 b/src/advance_field_mod.F90
index 7605db7a40e7a4441ce71d1e4cd79c96b3c4e68f..fd203489d74f253ca510c3f643bd69ecf5746167 100644
--- a/src/advance_field_mod.F90
+++ b/src/advance_field_mod.F90
@@ -16,59 +16,63 @@ CONTAINS
SUBROUTINE advance_moments
- USE basic
- USE time_integration
- USE grid
- use prec_const
+ USE basic, ONLY: t0_adv_field,t1_adv_field,tc_adv_field, dt
+ USE grid, ONLY:local_na,local_np,local_nj,local_nky,local_nkx,local_nz,&
+ ngp, ngj, ngz, dmax, parray, jarray, dmax
USE model, ONLY: CLOS
use fields, ONLY: moments
use array, ONLY: moments_rhs
+ USE time_integration, ONLY: updatetlevel, A_E, b_E, ntimelevel
IMPLICIT NONE
- INTEGER :: p_int, j_int, ia, ip, ij
-
+ INTEGER :: ia, ip, ij, ikx,iky,iz, istage, ipi, iji, izi
CALL cpu_time(t0_adv_field)
- DO ia=ias,iae
- DO ip=ips,ipe
- p_int = parray(ip)
- DO ij=ijs,ije
- j_int = jarray(ij)
- IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmax))&
- CALL advance_field(moments(ia,ip,ij,ikys:ikye,ikxs:ikxe,izs:ize,:), moments_rhs(ia,ip,ij,ikys:ikye,ikxs:ikxe,izs:ize,:))
- ENDDO
- ENDDO
- ENDDO
- ! Execution time end
- CALL cpu_time(t1_adv_field)
- tc_adv_field = tc_adv_field + (t1_adv_field - t0_adv_field)
- END SUBROUTINE advance_moments
-
-
- SUBROUTINE advance_field( f, f_rhs )
-
- USE basic
- USE time_integration
- USE array
- USE grid
- use prec_const
- IMPLICIT NONE
-
- COMPLEX(dp), DIMENSION ( ikys:ikye, ikxs:ikxe, izs:ize, ntimelevel ) :: f
- COMPLEX(dp), DIMENSION ( ikys:ikye, ikxs:ikxe, izs:ize, ntimelevel ) :: f_rhs
- INTEGER :: istage
-
SELECT CASE (updatetlevel)
CASE(1)
- DO istage=1,ntimelevel
- f(ikys:ikye,ikxs:ikxe,izs:ize,1) = f(ikys:ikye,ikxs:ikxe,izs:ize,1) &
- + dt*b_E(istage)*f_rhs(ikys:ikye,ikxs:ikxe,izs:ize,istage)
- END DO
+ DO istage=1,ntimelevel
+ DO iz =1,local_nz
+ izi = iz+ngz/2
+ DO ikx =1,local_nkx
+ DO iky =1,local_nky
+ DO ij =1,local_nj
+ iji = ij+ngj/2
+ DO ip =1,local_np
+ ipi = ip+ngp/2
+ DO ia =1,local_na
+ IF((CLOS .NE. 1) .OR. (parray(ip)+2*jarray(ij) .LE. dmax))&
+ moments(ia,ipi,iji,iky,ikx,izi,1) = moments(ia,ipi,iji,iky,ikx,izi,1) &
+ + dt*b_E(istage)*moments_rhs(ia,ip,ij,iky,ikx,iz,istage)
+ END DO
+ END DO
+ END DO
+ END DO
+ END DO
+ END DO
+ END DO
CASE DEFAULT
- f(ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel) = f(ikys:ikye,ikxs:ikxe,izs:ize,1);
- DO istage=1,updatetlevel-1
- f(ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel) = f(ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel) + &
- dt*A_E(updatetlevel,istage)*f_rhs(ikys:ikye,ikxs:ikxe,izs:ize,istage)
- END DO
+ moments(:,:,:,:,:,:,updatetlevel) = moments(:,:,:,:,:,:,1);
+ DO istage=1,ntimelevel-1
+ DO iz =1,local_nz
+ izi = iz+ngz/2
+ DO ikx =1,local_nkx
+ DO iky =1,local_nky
+ DO ij =1,local_nj
+ iji = ij+ngj/2
+ DO ip =1,local_np
+ ipi = ip+ngp/2
+ DO ia =1,local_na
+ IF((CLOS .NE. 1) .OR. (parray(ip)+2*jarray(ij) .LE. dmax))&
+ moments(ia,ipi,iji,iky,ikx,izi,updatetlevel) = moments(ia,ipi,iji,iky,ikx,izi,updatetlevel) + &
+ dt*A_E(updatetlevel,istage)*moments_rhs(ia,ip,ij,iky,ikx,iz,istage)
+ END DO
+ END DO
+ END DO
+ END DO
+ END DO
+ END DO
+ END DO
END SELECT
- END SUBROUTINE advance_field
-
+ ! Execution time end
+ CALL cpu_time(t1_adv_field)
+ tc_adv_field = tc_adv_field + (t1_adv_field - t0_adv_field)
+ END SUBROUTINE advance_moments
END MODULE advance_field_routine
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 7b46d23fbf72ba39cf1dbf361ee6c1e77f0c3351..d6f1fdeb5dadec560e7a5c95c85bc6a91c45924c 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -46,7 +46,7 @@ MODULE array
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_poisson_op
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_ampere_op
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_pol_ion
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: HF_phi_correction_operator
+ ! REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: HF_phi_correction_operator
! Gyrocenter density for electron and ions (ia,iky,ikx,iz)
COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Na00
diff --git a/src/auxval.F90 b/src/auxval.F90
index 3b2e27eebb1e88e8706084904b5e8571db11ecc1..d044790125f00c8b361a01e09452db2f7a3347b3 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -17,7 +17,7 @@ subroutine auxval
IF (my_id .EQ. 0) WRITE(*,*) '=== Set auxiliary values ==='
! Init the grids
- CALL set_grids(shear,Npol) ! radial modes (MPI distributed by FFTW)
+ CALL set_grids(shear,Npol,LINEARITY,N_HD,EM,Na) ! radial modes (MPI distributed by FFTW)
CALL memory ! Allocate memory for global arrays
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 9013db29c6d1b28e95481eaf8e633d422a8f30db..e8953b546d24df99f32e0f6472540ea4cc398742 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -142,13 +142,13 @@ CONTAINS
SUBROUTINE find_input_file
USE parallel, ONLY: my_id
IMPLICIT NONE
- CHARACTER(len=32) :: str, input_file
+ CHARACTER(len=32) :: str_, input_file
INTEGER :: nargs, fileid, l, ierr
LOGICAL :: mlexist
nargs = COMMAND_ARGUMENT_COUNT()
IF((nargs .EQ. 1) .OR. (nargs .EQ. 4)) THEN
- CALL GET_COMMAND_ARGUMENT(nargs, str, l, ierr)
- READ(str(1:l),'(i3)') fileid
+ CALL GET_COMMAND_ARGUMENT(nargs, str_, l, ierr)
+ READ(str_(1:l),'(i3)') fileid
WRITE(input_file,'(a,a1,i2.2,a3)') 'fort','_',fileid,'.90'
INQUIRE(file=input_file, exist=mlexist)
@@ -215,20 +215,20 @@ CONTAINS
END SUBROUTINE display_h_min_s
!================================================================================
- function str_dp(k) result( str )
+ function str_dp(k) result( str_ )
! "Convert an integer to string."
REAL(dp), intent(in) :: k
- character(len=20):: str
- write (str, *) k
- str = adjustl(str)
+ character(len=10):: str_
+ write (str_, "(G10.2)") k
+ str_ = adjustl(str_)
end function str_dp
- function str_int(k) result( str )
+ function str_int(k) result( str_ )
! "Convert an integer to string."
integer, intent(in) :: k
- character(len=20) :: str
- write (str, *) k
- str = adjustl(str)
+ character(len=10) :: str_
+ write (str_, "(i2.2)") k
+ str_ = adjustl(str_)
end function str_int
! To allocate arrays of doubles, integers, etc. at run time
diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
index a3ce639a86e0ea2b8f404c08a9b8bbd93805f2a7..d30b3a8c5299564a0a9ff20b8b29a8e514a4d3f1 100644
--- a/src/calculus_mod.F90
+++ b/src/calculus_mod.F90
@@ -38,7 +38,7 @@ SUBROUTINE grad_z(target,local_Nz,Ngz,inv_deltaz,f,ddzf)
CASE(2)
CALL grad_z_e2o(local_Nz,Ngz,inv_deltaz,f,ddzf)
CASE DEFAULT ! No staggered grid -> usual centered finite differences
- DO iz = Ngz/2+1,Ngz/2+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)
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index b68c5c031fc82d4d9753c3625a17ce6962b5931a..5c1585e05c3275dc4d2f36997887d4cfa8b6c978 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -112,19 +112,19 @@ SUBROUTINE update_ghosts_z_mom
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
!!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!!
! Send the last local moment to fill the -1 neighbour ghost
- DO ig=1,ngz
+ DO ig=1,ngz/2
CALL mpi_sendrecv(moments(:,:,:,:,:,last-(ig-1),updatetlevel),count,MPI_DOUBLE_COMPLEX,nbr_U,24+ig, & ! Send to Up the last
buff_pjxy_zBC(:,:,:,:,:,-ig),count,MPI_DOUBLE_COMPLEX,nbr_D,24+ig, & ! Recieve from Down the first-1
comm0, status, ierr)
ENDDO
!!!!!!!!!!! Send ghost to down neighbour !!!!!!!!!!!!!!!!!!!!!!
- DO ig=1,ngz
+ DO ig=1,ngz/2
CALL mpi_sendrecv(moments(:,:,:,:,:,first+(ig-1),updatetlevel),count,MPI_DOUBLE_COMPLEX,nbr_D,26+ig, & ! Send to Up the last
buff_pjxy_zBC(:,:,:,:,:,ig),count,MPI_DOUBLE_COMPLEX,nbr_U,26+ig, & ! Recieve from Down the first-1
comm0, status, ierr)
ENDDO
ELSE !No parallel (just copy)
- DO ig=1,ngz
+ DO ig=1,ngz/2
buff_pjxy_zBC(:,:,:,:,:,-ig) = moments(:,:,:,:,:, last-(ig-1),updatetlevel)
buff_pjxy_zBC(:,:,:,:,:, ig) = moments(:,:,:,:,:,first+(ig-1),updatetlevel)
ENDDO
@@ -135,11 +135,11 @@ SUBROUTINE update_ghosts_z_mom
! Exchanging the modes that have a periodic pair (from sheared BC)
IF (ikxBC_L .NE. -99) THEN
! Fill the lower ghosts cells with the buffer value (upper cells from LEFT process)
- DO ig=1,ngz
+ DO ig=1,ngz/2
moments(:,:,:,iky,ikx,first-ig,updatetlevel) = pb_phase_L(iky)*buff_pjxy_zBC(:,:,:,iky,ikxBC_L,-ig)
ENDDO
ELSE
- DO ig=1,ngz
+ DO ig=1,ngz/2
moments(:,:,:,iky,ikx,first-ig,updatetlevel) = 0._dp
ENDDO
ENDIF
@@ -147,11 +147,11 @@ SUBROUTINE update_ghosts_z_mom
! Exchanging the modes that have a periodic pair (from sheared BC)
IF (ikxBC_R .NE. -99) THEN
! Fill the upper ghosts cells with the buffer value (lower cells from RIGHT process)
- DO ig=1,ngz
+ DO ig=1,ngz/2
moments(:,:,:,iky,ikx,last+ig,updatetlevel) = pb_phase_R(iky)*buff_pjxy_zBC(:,:,:,iky,ikxBC_R,ig)
ENDDO
ELSE
- DO ig=1,ngz
+ DO ig=1,ngz/2
moments(:,:,:,iky,ikx,last+ig,updatetlevel) = 0._dp
ENDDO
ENDIF
@@ -174,7 +174,7 @@ SUBROUTINE update_ghosts_z_3D(field)
IF (num_procs_z .GT. 1) THEN
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
!!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!!
- DO ig = 1,ngz
+ DO ig = 1,ngz/2
CALL mpi_sendrecv( field(:,:,last-(ig-1)), count, MPI_DOUBLE_COMPLEX, nbr_U, 30+ig, & ! Send to Up the last
buff_xy_zBC(:,:,-ig), count, MPI_DOUBLE_COMPLEX, nbr_D, 30+ig, & ! Receive from Down the first-1
comm0, status, ierr)
@@ -187,7 +187,7 @@ SUBROUTINE update_ghosts_z_3D(field)
ENDDO
ELSE
! no parallelization so just copy last cell into first ghosts and vice versa
- DO ig = 1,ngz
+ DO ig = 1,ngz/2
buff_xy_zBC(:,:,-ig) = field(:,:,last-(ig-1))
buff_xy_zBC(:,:, ig) = field(:,:,first+(ig-1))
ENDDO
@@ -196,22 +196,22 @@ SUBROUTINE update_ghosts_z_3D(field)
DO ikx = 1,local_nkx
ikxBC_L = ikx_zBC_L(iky,ikx);
IF (ikxBC_L .GT. 0) THEN ! Exchanging the modes that have a periodic pair (a)
- DO ig = 1,ngz
+ DO ig = 1,ngz/2
field(iky,ikx,first-ig) = pb_phase_L(iky)*buff_xy_zBC(iky,ikxBC_L,-ig)
ENDDO
ELSE
- DO ig = 1,ngz
+ DO ig = 1,ngz/2
field(iky,ikx,first-ig) = 0._dp
ENDDO
ENDIF
ikxBC_R = ikx_zBC_R(iky,ikx);
IF (ikxBC_R .GT. 0) THEN ! Exchanging the modes that have a periodic pair (a)
! last+1 gets first
- DO ig = 1,ngz
+ DO ig = 1,ngz/2
field(iky,ikx,last+ig) = pb_phase_R(iky)*buff_xy_zBC(iky,ikxBC_R,ig)
ENDDO
ELSE
- DO ig = 1,ngz
+ DO ig = 1,ngz/2
field(iky,ikx,last+ig) = 0._dp
ENDDO
ENDIF
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index ed59684e427abcf0e1fae3839cfe4379423a6b08..a344e07923e744178d6bf7412218949d4b638a73 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -48,14 +48,14 @@ MODULE grid
INTEGER, PUBLIC, PROTECTED :: total_nkx
INTEGER, PUBLIC, PROTECTED :: total_nz
! Local numbers of points (without ghosts)
- INTEGER, PUBLIC, PROTECTED :: local_Na
- INTEGER, PUBLIC, PROTECTED :: local_Np
- INTEGER, PUBLIC, PROTECTED :: local_Nj
- INTEGER, PUBLIC, PROTECTED :: local_Nky
- INTEGER, PUBLIC, PROTECTED :: local_Nkx
- INTEGER, PUBLIC, PROTECTED :: local_Nz
- INTEGER, PUBLIC, PROTECTED :: local_Nkp
- INTEGER, PUBLIC, PROTECTED :: Ngp, Ngj, Ngx, Ngy, Ngz ! number of ghosts points
+ INTEGER, PUBLIC, PROTECTED :: local_na
+ INTEGER, PUBLIC, PROTECTED :: local_np
+ INTEGER, PUBLIC, PROTECTED :: local_nj
+ INTEGER, PUBLIC, PROTECTED :: local_nky
+ INTEGER, PUBLIC, PROTECTED :: local_nkx
+ INTEGER, PUBLIC, PROTECTED :: local_nz
+ INTEGER, PUBLIC, PROTECTED :: local_nkp
+ INTEGER, PUBLIC, PROTECTED :: ngp, ngj, ngx, ngy, ngz ! number of ghosts points
INTEGER, PUBLIC, PROTECTED :: Nzgrid ! one or two depending on the staggered grid option
! Local offsets
INTEGER, PUBLIC, PROTECTED :: local_na_offset
@@ -70,8 +70,8 @@ MODULE grid
! Grid spacing and limits
REAL(dp), PUBLIC, PROTECTED :: deltap, pp2, deltaz, inv_deltaz
REAL(dp), PUBLIC, PROTECTED :: deltakx, deltaky, kx_max, ky_max, kx_min, ky_min!, kp_max
- REAL(dp), PUBLIC, PROTECTED :: local_pmin, local_pmax
- REAL(dp), PUBLIC, PROTECTED :: local_jmin, local_jmax
+ INTEGER , PUBLIC, PROTECTED :: local_pmin, local_pmax
+ INTEGER , PUBLIC, PROTECTED :: local_jmin, local_jmax
REAL(dp), PUBLIC, PROTECTED :: local_kymin, local_kymax
REAL(dp), PUBLIC, PROTECTED :: local_kxmin, local_kxmax
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: local_zmin, local_zmax
@@ -141,13 +141,17 @@ CONTAINS
! The other grids are simply
! |_|_|_|_|
! 1 2 3 4
- SUBROUTINE set_grids(shear,Npol)
- USE model, ONLY: LINEARITY
+ SUBROUTINE set_grids(shear,Npol,LINEARITY,N_HD,EM,Na)
USE fourier, ONLY: init_grid_distr_and_plans
REAL(dp), INTENT(IN) :: shear
INTEGER, INTENT(IN) :: Npol
- CALL set_agrid
- CALL set_pgrid
+ CHARACTER(len=*), INTENT(IN) :: LINEARITY
+ INTEGER, INTENT(IN) :: N_HD
+ LOGICAL, INTENT(IN) :: EM
+ INTEGER, INTENT(IN) :: Na
+ CALL set_agrid(Na)
+ CALL set_pgrid(EM)
+ CALL set_jgrid
!! Parallel distribution of kx ky grid
IF (LINEARITY .NE. 'linear') THEN
IF (my_id .EQ. 0) write(*,*) 'FFTW3 y-grid distribution'
@@ -156,8 +160,8 @@ CONTAINS
CALL init_1Dgrid_distr
IF (my_id .EQ. 0) write(*,*) 'Manual y-grid distribution'
ENDIF
- CALL set_kygrid
- CALL set_kxgrid(shear,Npol)
+ CALL set_kygrid(LINEARITY,N_HD)
+ CALL set_kxgrid(shear,Npol,LINEARITY,N_HD)
CALL set_zgrid (Npol)
END SUBROUTINE set_grids
@@ -170,9 +174,9 @@ CONTAINS
if (rank_ky .EQ. num_procs_ky-1) local_nky_ptr = (Ny/2+1)-local_nky_ptr_offset
END SUBROUTINE init_1Dgrid_distr
- SUBROUTINE set_agrid ! you're a sorcerer harry
- USE model, ONLY: Na
+ SUBROUTINE set_agrid(Na) ! you're a sorcerer harry
IMPLICIT NONE
+ INTEGER, INTENT(IN) :: Na
ias = 1
iae = Na
total_Na = Na
@@ -180,11 +184,11 @@ CONTAINS
local_Na_offset = ias - 1
END SUBROUTINE
- SUBROUTINE set_pgrid
+ SUBROUTINE set_pgrid(EM)
USE prec_const
- USE model, ONLY: beta ! To know if we solve ampere or not and put odd p moments
USE parallel, ONLY: num_procs_p, rank_p
IMPLICIT NONE
+ LOGICAL, INTENT(IN) :: EM
INTEGER :: ip, ig
! If no parallel dim (Nz=1) and no EM effects (beta=0), the moment hierarchy
!! is separable between odds and even P and since the energy is injected in
@@ -192,7 +196,7 @@ CONTAINS
!! simulating the odd p which will only be damped.
!! We define in this case a grid Parray = 0,2,4,...,Pmax i.e. deltap = 2
!! instead of 1 to spare computation
- IF((Nz .EQ. 1) .AND. (beta .EQ. 0._dp)) THEN
+ IF((Nz .EQ. 1) .AND. .NOT. EM) THEN
deltap = 2
Ngp = 2 ! two ghosts cells for p+/-2 only
pp2 = 1 ! index p+2 is ip+1
@@ -246,7 +250,7 @@ CONTAINS
pmax_dp = real(pmax,dp)
diff_p_coeff = pmax_dp*(1._dp/pmax_dp)**6
! Overwrite SOLVE_AMPERE flag if beta is zero
- IF(beta .EQ. 0._dp) THEN
+ IF(.NOT. EM) THEN
SOLVE_AMPERE = .FALSE.
ENDIF
END SUBROUTINE set_pgrid
@@ -262,7 +266,7 @@ CONTAINS
! Build the full grids on process 0 to diagnose it without comm
ALLOCATE(jarray_full(total_nj))
! J
- DO ij = 1,jmax+1; jarray_full(ij) = (ij-1); END DO
+ DO ij = 1,total_nj; jarray_full(ij) = (ij-1); END DO
! Indices of local data
ijs = 1; ije = jmax + 1
! Local number of J
@@ -272,8 +276,8 @@ CONTAINS
DO ij = 1+ngj/2,local_nj+ngj/2
jarray(ij) = ij-1-ngj/2+local_nj_offset
END DO
- local_jmax = jarray(local_np+ngp/2)
- local_jmin = jarray(1+ngp/2)
+ local_jmax = jarray(local_np+ngj/2)
+ local_jmin = jarray(1+ngj/2)
! Fill the ghosts
DO ig = 1,ngj/2
jarray(ig) = local_jmin-ngj/2+(ig-1)
@@ -288,10 +292,11 @@ CONTAINS
END DO
END SUBROUTINE set_jgrid
- SUBROUTINE set_kygrid
+ SUBROUTINE set_kygrid(LINEARITY,N_HD)
USE prec_const
- USE model, ONLY: LINEARITY, N_HD
IMPLICIT NONE
+ CHARACTER(len=*), INTENT(IN) ::LINEARITY
+ INTEGER, INTENT(IN) :: N_HD
INTEGER :: iky, ikyo
Nky = Ny/2+1 ! Defined only on positive kx since fields are real
! Grid spacings
@@ -362,12 +367,13 @@ CONTAINS
ENDIF
END SUBROUTINE set_kygrid
- SUBROUTINE set_kxgrid(shear,Npol)
+ SUBROUTINE set_kxgrid(shear,Npol,LINEARITY,N_HD)
USE prec_const
- USE model, ONLY: LINEARITY, N_HD
IMPLICIT NONE
REAL(dp), INTENT(IN) :: shear
INTEGER, INTENT(IN) :: Npol
+ CHARACTER(len=*), INTENT(IN) ::LINEARITY
+ INTEGER, INTENT(IN) :: N_HD
INTEGER :: ikx, ikxo
REAL :: Lx_adapted
IF(shear .GT. 0) THEN
@@ -489,7 +495,6 @@ CONTAINS
SUBROUTINE set_zgrid(Npol)
USE prec_const
- USE model, ONLY: mu_z
USE parallel, ONLY: num_procs_z, rank_z
IMPLICIT NONE
REAL :: grid_shift, Lz, zmax, zmin
@@ -501,11 +506,7 @@ CONTAINS
deltaz = Lz/REAL(Nz,dp)
inv_deltaz = 1._dp/deltaz
! Parallel hyperdiffusion coefficient
- IF(mu_z .GT. 0) THEN
- diff_dz_coeff = (deltaz/2._dp)**4 ! adaptive fourth derivative (~GENE)
- ELSE
- diff_dz_coeff = -1._dp ! non adaptive (negative sign to compensate mu_z neg)
- ENDIF
+ diff_dz_coeff = (deltaz/2._dp)**4 ! adaptive fourth derivative (~GENE)
IF (SG) THEN
CALL speak('--2 staggered z grids--')
grid_shift = deltaz/2._dp
diff --git a/src/memory.F90 b/src/memory.F90
index 4263680f9a6f45fe75a5759d88d00f0badebc971..c9d856d17fb19c3528814c32d2de6fc76acceb75 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -17,7 +17,7 @@ SUBROUTINE memory
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(HF_phi_correction_operator, 1,local_nky, 1,local_nkx, 1,local_Nz)
!Moments related arrays
CALL allocate_array( Na00, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz)
@@ -32,9 +32,9 @@ SUBROUTINE memory
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+Ngz)
- 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+Ngz)
- 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+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)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 599fe2c4afa771363b7dee4f6795e82f702ba6de..db3a6c85cd59fc1df3859ae77ef922ade60777f9 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -51,10 +51,10 @@ SUBROUTINE compute_moments_eq_rhs
p_int = parray(ipi) ! Hermite degree
eo = min(nzgrid,MODULO(p_int,2)+1) ! Indicates if we are on odd or even z grid
kperp2= kparray(iky,ikx,izi,eo)**2
- Napj = moments(ia,ipi,iji,iky,ikx,izi,updatetlevel)
RHS = 0._dp
! Species loop
a:DO ia = 1,local_na
+ Napj = moments(ia,ipi,iji,iky,ikx,izi,updatetlevel)
IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmax)) THEN ! for the closure scheme
!! Compute moments_ mixing terms
! Perpendicular dynamic
@@ -72,17 +72,17 @@ SUBROUTINE compute_moments_eq_rhs
Mperp = imagu*Ckxky(iky,ikx,iz,eo)*(Tnapj + Tnapp2j + Tnapm2j + Tnapjp1 + Tnapjm1)
! Parallel dynamic
! ddz derivative for Landau damping term
- Ldamp = xnapp1j(ia,ip) * ddz_napj(ia,ipi+1,ij,iky,ikx,izi) &
- + xnapm1j(ia,ip) * ddz_napj(ia,ipi-1,ij,iky,ikx,izi)
+ Ldamp = xnapp1j(ia,ip) * ddz_napj(ia,ipi+1,iji,iky,ikx,iz) &
+ + xnapm1j(ia,ip) * ddz_napj(ia,ipi-1,iji,iky,ikx,iz)
! Mirror terms
- Tnapp1j = ynapp1j (ia,ip,ij) * interp_napj(ia,ipi+1,ij ,iky,ikx,izi)
- Tnapp1jm1 = ynapp1jm1(ia,ip,ij) * interp_napj(ia,ipi+1,ij-1,iky,ikx,izi)
- Tnapm1j = ynapm1j (ia,ip,ij) * interp_napj(ia,ipi-1,ij ,iky,ikx,izi)
- Tnapm1jm1 = ynapm1jm1(ia,ip,ij) * interp_napj(ia,ipi-1,ij-1,iky,ikx,izi)
+ Tnapp1j = ynapp1j (ia,ip,ij) * interp_napj(ia,ipi+1,iji ,iky,ikx,iz)
+ Tnapp1jm1 = ynapp1jm1(ia,ip,ij) * interp_napj(ia,ipi+1,iji-1,iky,ikx,iz)
+ Tnapm1j = ynapm1j (ia,ip,ij) * interp_napj(ia,ipi-1,iji ,iky,ikx,iz)
+ Tnapm1jm1 = ynapm1jm1(ia,ip,ij) * interp_napj(ia,ipi-1,iji-1,iky,ikx,iz)
! Trapping terms
- Unapm1j = znapm1j (ia,ip,ij) * interp_napj(ia,ipi-1,ij ,iky,ikx,izi)
- Unapm1jp1 = znapm1jp1(ia,ip,ij) * interp_napj(ia,ipi-1,ij+1,iky,ikx,izi)
- Unapm1jm1 = znapm1jm1(ia,ip,ij) * interp_napj(ia,ipi-1,ij-1,iky,ikx,izi)
+ Unapm1j = znapm1j (ia,ip,ij) * interp_napj(ia,ipi-1,iji ,iky,ikx,iz)
+ Unapm1jp1 = znapm1jp1(ia,ip,ij) * interp_napj(ia,ipi-1,iji+1,iky,ikx,iz)
+ Unapm1jm1 = znapm1jm1(ia,ip,ij) * interp_napj(ia,ipi-1,iji-1,iky,ikx,iz)
! sum the parallel forces
Fmir = dlnBdz(iz,eo)*(Tnapp1j + Tnapp1jm1 + Tnapm1j + Tnapm1jm1 +&
Unapm1j + Unapm1jp1 + Unapm1jm1)
@@ -124,7 +124,7 @@ SUBROUTINE compute_moments_eq_rhs
-mu_x*diff_kx_coeff*kx**N_HD*Napj &
-mu_y*diff_ky_coeff*ky**N_HD*Napj &
! Numerical parallel hyperdiffusion "mu_z*ddz**4" see Pueschel 2010 (eq 25)
- -mu_z*diff_dz_coeff*ddzND_napj(ia,ipi,iji,iky,ikx,izi)
+ -mu_z*diff_dz_coeff*ddzND_napj(ia,ipi,iji,iky,ikx,iz)
!! Velocity space dissipation (should be implemented somewhere else)
SELECT CASE(HYP_V)
CASE('hypcoll') ! GX like Hermite hypercollisions see Mandell et al. 2023 (eq 3.23), unadvised to use it
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index e62df88068e9092051e6352dce76e41cdad3702c..189ac0a9d41bb328753275524243cb1c9d8d9852 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -71,7 +71,7 @@ END SUBROUTINE build_dv4Hp_table
!******************************************************************************!
SUBROUTINE evaluate_kernels
USE basic
- USE array, ONLY : kernel, HF_phi_correction_operator
+ 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 species, ONLY : sigma2_tau_o2
USE prec_const, ONLY: dp
@@ -101,20 +101,20 @@ DO ia = 1,local_Na
ENDDO
ENDDO
ENDDO
- !! Correction term for the evaluation of the heat flux
- HF_phi_correction_operator(:,:,:) = &
- 2._dp * Kernel(ia,1,:,:,:,1) &
- -1._dp * Kernel(ia,2,:,:,:,1)
-
- DO ij = 1, local_Nj
- j_int = jarray(ij)
- j_dp = REAL(j_int,dp)
- HF_phi_correction_operator(:,:,:) = HF_phi_correction_operator(:,:,:) &
- - Kernel(ia,ij,:,:,:,1) * (&
- 2._dp*(j_dp+1.5_dp) * Kernel(ia,ij ,:,:,:,1) &
- - (j_dp+1.0_dp) * Kernel(ia,ij+1,:,:,:,1) &
- - j_dp * Kernel(ia,ij-1,:,:,:,1))
- ENDDO
+ ! !! Correction term for the evaluation of the heat flux
+ ! HF_phi_correction_operator(:,:,:) = &
+ ! 2._dp * Kernel(ia,1,:,:,:,1) &
+ ! -1._dp * Kernel(ia,2,:,:,:,1)
+ !
+ ! DO ij = 1, local_Nj
+ ! j_int = jarray(ij)
+ ! j_dp = REAL(j_int,dp)
+ ! HF_phi_correction_operator(:,:,:) = HF_phi_correction_operator(:,:,:) &
+ ! - Kernel(ia,ij,:,:,:,1) * (&
+ ! 2._dp*(j_dp+1.5_dp) * Kernel(ia,ij ,:,:,:,1) &
+ ! - (j_dp+1.0_dp) * Kernel(ia,ij+1,:,:,:,1) &
+ ! - j_dp * Kernel(ia,ij-1,:,:,:,1))
+ ! ENDDO
ENDDO
END SUBROUTINE evaluate_kernels
@@ -134,7 +134,7 @@ SUBROUTINE evaluate_poisson_op
USE basic
USE array, ONLY : kernel, inv_poisson_op, inv_pol_ion
USE grid, ONLY : local_Na, local_nkx, local_nky, local_nz,&
- kxarray, kyarray, jmax
+ kxarray, kyarray, local_nj, ngj, ngz, ieven
USE species, ONLY : q2_tau
USE model, ONLY : ADIAB_E, tau_e
USE prec_const, ONLY: dp
@@ -155,12 +155,13 @@ SUBROUTINE evaluate_poisson_op
pol_tot = 0._dp ! total polarisation term
pol_ion = 0._dp ! sum of ion polarisation term
! loop over n only up to the max polynomial degree
- DO in=1,jmax+1
- pol_tot = pol_tot + q2_tau(ia)*kernel(ia,in,iky,ikx,iz,0)**2 ! ... sum recursively ...
- pol_ion = pol_ion + q2_tau(ia)*kernel(ia,in,iky,ikx,iz,0)**2 !
+ DO in=1,local_nj
+ pol_tot = pol_tot + q2_tau(ia)*kernel(ia,in+ngj/2,iky,ikx,iz+ngz/2,ieven)**2 ! ... sum recursively ...
+ pol_ion = pol_ion + q2_tau(ia)*kernel(ia,in+ngj/2,iky,ikx,iz+ngz/2,ieven)**2 !
END DO
operator = operator + q2_tau(ia) - pol_tot
ENDDO
+ operator_ion = operator
IF(ADIAB_E) THEN ! Adiabatic model
pol_tot = pol_tot + 1._dp/tau_e - 1._dp
ENDIF
@@ -227,17 +228,17 @@ SUBROUTINE compute_lin_coeff
USE species, ONLY: k_T, k_N, tau, q, sqrtTau_q, tau_q
USE model, ONLY: k_cB, k_gB
USE prec_const, ONLY: dp, SQRT2, SQRT3
- USE grid, ONLY: parray, jarray, local_na, local_np, local_nj
+ USE grid, ONLY: parray, jarray, local_na, local_np, local_nj, ngj, ngp
INTEGER :: ia,ip,ij,p_int, j_int ! polynom. dagrees
REAL(dp) :: p_dp, j_dp
!! linear coefficients for moment RHS !!!!!!!!!!
DO ia = 1,local_na
DO ip = 1, local_np
- p_int= parray(ip) ! Hermite degree
+ p_int= parray(ip+ngp/2) ! Hermite degree
p_dp = REAL(p_int,dp) ! REAL of Hermite degree
DO ij = 1, local_nj
- j_int= jarray(ij) ! Laguerre degree
+ j_int= jarray(ij+ngj/2) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
! All Napj terms
xnapj(ia,ip,ij) = tau(ia)/q(ia)*(k_cB*(2._dp*p_dp + 1._dp) &
@@ -254,7 +255,7 @@ SUBROUTINE compute_lin_coeff
ENDDO
ENDDO
DO ip = 1, local_np
- p_int= parray(ip) ! Hermite degree
+ p_int= parray(ip+ngp/2) ! Hermite degree
p_dp = REAL(p_int,dp) ! REAL of Hermite degree
! Landau damping coefficients (ddz napj term)
xnapp1j(ia,ip) = sqrtTau_q(ia) * SQRT(p_dp+1._dp)
@@ -264,7 +265,7 @@ SUBROUTINE compute_lin_coeff
xnapm2j(ia,ip) = tau_q(ia) * k_cB * SQRT( p_dp *(p_dp - 1._dp))
ENDDO
DO ij = 1, local_nj
- j_int= jarray(ij) ! Laguerre degree
+ j_int= jarray(ij+ngj/2) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
! Magnetic gradient term
xnapjp1(ia,ij) = -tau_q(ia) * k_gB * (j_dp + 1._dp)
@@ -273,9 +274,9 @@ SUBROUTINE compute_lin_coeff
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! ES linear coefficients for moment RHS !!!!!!!!!!
DO ip = 1, local_np
- p_int= parray(ip) ! Hermite degree
+ p_int= parray(ip+ngp/2) ! Hermite degree
DO ij = 1, local_nj
- j_int= jarray(ij) ! REALof Laguerre degree
+ j_int= jarray(ij+ngj/2) ! REALof Laguerre degree
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 0) THEN ! kronecker p0
@@ -294,9 +295,9 @@ SUBROUTINE compute_lin_coeff
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Electromagnatic linear coefficients for moment RHS !!!!!!!!!!
DO ip = 1, local_np
- p_int= parray(ip) ! Hermite degree
+ p_int= parray(ip+ngp/2) ! Hermite degree
DO ij = 1, local_nj
- j_int= jarray(ij) ! REALof Laguerre degree
+ j_int= jarray(ij+ngj/2) ! REALof Laguerre degree
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
IF (p_int .EQ. 1) THEN ! kronecker p1
xpsij (ia,ip,ij) = +(k_N(ia) + (2._dp*j_dp+1._dp)*k_T(ia))* sqrtTau_q(ia)
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 49aaf22ee6b76f1d65a06002077c4737f263e847..ecfe01c7bf0d873aa38be6c1b76d81e6555ce813 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -2,7 +2,7 @@ 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,&
+ parray,pmax,ip0, iodd, ieven,&
CONTAINSp0,ip1,CONTAINSp1,ip2,CONTAINSp2,ip3,CONTAINSp3,&
jarray,jmax,ij0, dmax,&
kyarray, AA_y,&
@@ -233,6 +233,8 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
!
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
@@ -281,15 +283,23 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
p_int = parray(ip)
eo = MODULO(p_int,2)+1 ! Indicates if we are on even or odd z grid
! Compute z first derivative
- CALL grad_z(eo,local_nz,ngz,inv_deltaz,nadiab_moments(ia,ip,ij,iky,ikx,:),ddz_napj(ia,ip,ij,iky,ikx,:))
+ 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
- CALL grad_z4(local_nz,ngz,inv_deltaz,nadiab_moments(ia,ip,ij,iky,ikx,:),ddzND_Napj(ia,ip,ij,iky,ikx,:))
+ 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
- CALL grad_z4(local_nz,ngz,inv_deltaz,moments(ia,ip,ij,iky,ikx,:,updatetlevel),ddzND_Napj(ia,ip,ij,iky,ikx,:))
+ 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
- CALL interp_z(eo,local_nz,ngz,nadiab_moments(ia,ip,ij,iky,ikx,1:local_nz+ngz), interp_napj(ia,ip,ij,iky,ikx,1:local_nz))
+ 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
@@ -369,7 +379,7 @@ SUBROUTINE compute_density
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,0) * moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
+ 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
@@ -391,7 +401,7 @@ SUBROUTINE compute_uperp
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,0) *&
+ 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_
@@ -414,7 +424,7 @@ SUBROUTINE compute_upar
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,1)*moments(ia,ip1,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
+ 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
@@ -440,7 +450,7 @@ SUBROUTINE compute_tperp
Tperp_ = 0._dp
DO ij = 1, local_nj
j_dp = REAL(ij-1,dp)
- Tperp_ = Tperp_ + kernel(ia,ij,iky,ikx,iz,0)*&
+ 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))
@@ -470,7 +480,7 @@ SUBROUTINE compute_Tpar
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,0)*&
+ 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
diff --git a/testcases/smallest_problem/fort.90 b/testcases/smallest_problem/fort.90
index 532f1db0ce64daf09fd240622434082c4c90c7c4..7c491e5dfd2df784fbb611b2bf1dda86084484a5 100644
--- a/testcases/smallest_problem/fort.90
+++ b/testcases/smallest_problem/fort.90
@@ -1,7 +1,7 @@
&BASIC
- nrun = 10
+ nrun = 500
dt = 0.01
- tmax = 1
+ tmax = 5
maxruntime = 356400
job2load = -1
/