diff --git a/Makefile b/Makefile
index 42d2305007a76e0c725e516c5f2a69f640dbb0bb..a0820e0a2c96c441c54c877160d970a59430491f 100644
--- a/Makefile
+++ b/Makefile
@@ -2,7 +2,7 @@ include local/dirs.inc
include local/make.inc
-EXEC = $(BINDIR)/monli1D
+EXEC = $(BINDIR)/helaz
F90 = mpif90
@@ -38,21 +38,21 @@ $(OBJDIR)/diagnose.o : src/srcinfo.h
FOBJ=$(OBJDIR)/advance_field.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o $(OBJDIR)/basic_mod.o \
$(OBJDIR)/control.o $(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o \
-$(OBJDIR)/fft_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/gradients.o $(OBJDIR)/inital.o \
+$(OBJDIR)/fields_mod.o $(OBJDIR)/inital.o \
$(OBJDIR)/initial_par_mod.o $(OBJDIR)/main.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o \
-$(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/numerics.o $(OBJDIR)/poisson.o $(OBJDIR)/ppexit.o \
-$(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/readinputs.o $(OBJDIR)/space_grid_mod.o \
-$(OBJDIR)/stepon.o $(OBJDIR)/temp_eq_rhs.o $(OBJDIR)/tesend.o $(OBJDIR)/theta_eq_rhs.o \
-$(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/vpar_eq_rhs.o
+$(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o $(OBJDIR)/ppexit.o \
+$(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/readinputs.o $(OBJDIR)/fourier_grid_mod.o \
+$(OBJDIR)/stepon.o $(OBJDIR)/tesend.o \
+$(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
$(EXEC): $(FOBJ)
$(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@
-$(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/space_grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
+$(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/advance_field.F90 -o $@
$(OBJDIR)/array_mod.o : src/array_mod.F90 $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/array_mod.F90 -o $@
-$(OBJDIR)/auxval.o : src/auxval.F90 $(OBJDIR)/fft_mod.o $(OBJDIR)/gradients.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o $(OBJDIR)/space_grid_mod.o
+$(OBJDIR)/auxval.o : src/auxval.F90 $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o $(OBJDIR)/fourier_grid_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/auxval.F90 -o $@
$(OBJDIR)/basic_mod.o : src/basic_mod.F90 $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/basic_mod.F90 -o $@
@@ -64,50 +64,38 @@ $(OBJDIR)/diagnostics_par_mod.o : src/diagnostics_par_mod.F90 $(OBJDIR)/prec_con
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/diagnostics_par_mod.F90 -o $@
$(OBJDIR)/endrun.o : src/endrun.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/endrun.F90 -o $@
-$(OBJDIR)/fft_mod.o : src/fft_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fft_mod.F90 -o $@
$(OBJDIR)/fields_mod.o : src/fields_mod.F90 $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fields_mod.F90 -o $@
-$(OBJDIR)/gradients.o : src/gradients.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/space_grid_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/gradients.F90 -o $@
-$(OBJDIR)/inital.o : src/inital.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/space_grid_mod.o $(OBJDIR)/time_integration_mod.o
+$(OBJDIR)/inital.o : src/inital.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/inital.F90 -o $@
$(OBJDIR)/initial_par_mod.o : src/initial_par_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/initial_par_mod.F90 -o $@
$(OBJDIR)/main.o : src/main.F90 $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/main.F90 -o $@
-$(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/space_grid_mod.o
+$(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/fourier_grid_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/memory.F90 -o $@
$(OBJDIR)/model_mod.o : src/model_mod.F90 $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/model_mod.F90 -o $@
-$(OBJDIR)/moments_eq_rhs.o : src/moments_eq_rhs.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/space_grid_mod.o $(OBJDIR)/time_integration_mod.o
+$(OBJDIR)/moments_eq_rhs.o : src/moments_eq_rhs.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs.F90 -o $@
-$(OBJDIR)/numerics.o : src/numerics.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/space_grid_mod.o $(OBJDIR)/time_integration_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/numerics.F90 -o $@
-$(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/space_grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
+$(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/poisson.F90 -o $@
$(OBJDIR)/ppexit.o : src/ppexit.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppexit.F90 -o $@
-$(OBJDIR)/ppinit.o : src/ppinit.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/space_grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
+$(OBJDIR)/ppinit.o : src/ppinit.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppinit.F90 -o $@
$(OBJDIR)/prec_const_mod.o : src/prec_const_mod.F90
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/prec_const_mod.F90 -o $@
-$(OBJDIR)/readinputs.o : src/readinputs.F90 $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/space_grid_mod.o $(OBJDIR)/time_integration_mod.o
+$(OBJDIR)/readinputs.o : src/readinputs.F90 $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/readinputs.F90 -o $@
-$(OBJDIR)/space_grid_mod.o : src/space_grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/space_grid_mod.F90 -o $@
-$(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/numerics.o $(OBJDIR)/space_grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o $(OBJDIR)/temp_eq_rhs.o $(OBJDIR)/theta_eq_rhs.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/vpar_eq_rhs.o $(OBJDIR)/model_mod.o
+$(OBJDIR)/fourier_grid_mod.o : src/fourier_grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_grid_mod.F90 -o $@
+$(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/stepon.F90 -o $@
-$(OBJDIR)/temp_eq_rhs.o : src/temp_eq_rhs.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/space_grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/temp_eq_rhs.F90 -o $@
$(OBJDIR)/tesend.o : src/tesend.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/tesend.F90 -o $@
-$(OBJDIR)/theta_eq_rhs.o : src/theta_eq_rhs.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/space_grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/theta_eq_rhs.F90 -o $@
$(OBJDIR)/time_integration_mod.o : src/time_integration_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/time_integration_mod.F90 -o $@
-$(OBJDIR)/utility_mod.o : src/utility_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/space_grid_mod.o
+$(OBJDIR)/utility_mod.o : src/utility_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
-$(OBJDIR)/vpar_eq_rhs.o : src/vpar_eq_rhs.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/space_grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/vpar_eq_rhs.F90 -o $@
diff --git a/src/advance_field.F90 b/src/advance_field.F90
index 9a968d1411b0d34177061690560dc7625c388d7f..0780a9315f55fb54378168ae60f673ab9fa762b7 100644
--- a/src/advance_field.F90
+++ b/src/advance_field.F90
@@ -23,12 +23,12 @@ CONTAINS
USE basic
USE time_integration
USE array
- USE space_grid
+ USE fourier_grid
use prec_const
IMPLICIT NONE
- real(dp), DIMENSION ( ikrs:ikre,ikzs:ikze,ntimelevel ) :: f
- real(dp), DIMENSION ( ikrs:ikre,ikzs:ikze,ntimelevel ) :: f_rhs
+ COMPLEX(dp), DIMENSION ( ikrs:ikre, ikzs:ikze, ntimelevel ) :: f
+ COMPLEX(dp), DIMENSION ( ikrs:ikre, ikzs:ikze, ntimelevel ) :: f_rhs
INTEGER :: ikr, ikz, istage
SELECT CASE (updatetlevel)
@@ -42,7 +42,7 @@ CONTAINS
END DO
CASE DEFAULT
DO ikz=ikzs,ikze
- f(ikr,ikz,updatetlevel) = f(ikz,1);
+ f(ikr,ikz,updatetlevel) = f(ikr,ikz,1);
DO ikr=ikrs,ikre
DO istage=1,updatetlevel-1
f(ikr,ikz,updatetlevel) = f(ikr,ikz,updatetlevel) + &
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index e0c0e655dbe77c6cb8e2ac96f88bcbdc84c9e98e..1fed4a4adba0903b215f27965941b8cf44a9aa31 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -4,64 +4,11 @@ MODULE array
use prec_const
implicit none
- ! Fields, on opposite grid
- real(dp), DIMENSION(:), ALLOCATABLE :: vpar_n
-
! Arrays to store the rhs, for time integration
- real(dp), DIMENSION(:,:,:), ALLOCATABLE :: theta_rhs, temp_rhs, vpar_rhs ! (ikr,ikz,updatetlevel)
- real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_rhs ! (ip,ij,ikr,ikz,updatetlevel)
- ! rhs, on opposite grid
- real(dp), DIMENSION(:,:), ALLOCATABLE :: temp_rhs_v, vpar_rhs_n !(ikr,ikz)
-
-
- real(dp), DIMENSION(:,:), ALLOCATABLE:: sqrt_exp_temp ! Often needed in 'rhs', compute it only once and store it
- real(dp), DIMENSION(:,:), ALLOCATABLE:: sqrt_exp_temp_v ! On opposite grid
-
- ! Spatial 1rst derivatives on respective grids
- ! real(dp), DIMENSION(:,:), ALLOCATABLE:: thetaz, sqrt_exp_tempz, vparz, phiz ! (ikr,ikz)
- ! real(dp), DIMENSION(:,:,:,:), ALLOCATABLE:: momentsz ! (ip,ij,ikr,ikz)
- ! Spatial 1rst derivatives on opposite grid
- ! real(dp), DIMENSION(:,:), ALLOCATABLE:: thetaz_v, sqrt_exp_tempz_v, vparz_n, phiz_v ! (ikr,ikz)
-
- ! ! For upwind scheme, other set of spatial derivatives, the "minus side"
- ! real(dp), DIMENSION(:), ALLOCATABLE:: thetazm, tempzm, vparzm, phizm
- ! real(dp), DIMENSION(:,:), ALLOCATABLE:: momentszm
-
- ! Spatial 2nd derivatives, for numerical diffusion, on respective grids
- ! real(dp), DIMENSION(:,:), ALLOCATABLE:: thetazz, tempzz, vparzz
- ! real(dp), DIMENSION(:,:,:), ALLOCATABLE:: momentszz
+ COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments_rhs ! (ipj,ikr,ikz,updatetlevel)
! Intermediate steps in rhs of equations
- real(dp), DIMENSION(:,:), ALLOCATABLE:: theta_Cpl, theta_Hpl
- real(dp), DIMENSION(:,:), ALLOCATABLE:: temp_Cpl, temp_Fpl, temp_Hpl, temp_Ipl
- real(dp), DIMENSION(:,:), ALLOCATABLE:: vpar_Cpl, vpar_Dpl, vpar_Epl, vpar_Fpl, vpar_Hpl
- real(dp), DIMENSION(:,:), ALLOCATABLE:: moments_Apl, moments_Bpl, moments_Cpl, moments_Dpl, moments_Epl, moments_Fpl, moments_Gpl, moments_Hpl
- ! moments_Apl = A^l
- ! moments_Bpl = sum_p [ B_p^l N_e^p ]
- ! moments_Cpl = sum_p [ C_p^l N_e^p ]
- ! moments_Dpl = sum_p [ D_p^l N_e^p ]
- ! moments_Epl = sum_p [ E_p^l N_e^p ]
- ! moments_Fpl = sum_p [ F_p^l N_e^p ]
- ! moments_Gpl = sum_p [ G_p^l N_e^p ]
- ! moments_Hpl = sum_p [ H_p^l N_e^p ]
- real(dp) :: moments_Ipl
- ! moments_Ipl = sum_p [ I_p^l d(N_e^p)/dz ]
-
-
- ! For poisson solver
- !TYPE(mumps_mat), SAVE :: laplace_smumps
- real(dp), DIMENSION(:,:), ALLOCATABLE :: laplace_rhs ! (ikr,ikz)
- real(dp), DIMENSION(:,:), ALLOCATABLE :: laplace_sol ! (ikr,ikz)
-
- ! For FFT
- double complex, DIMENSION(:,:), ALLOCATABLE :: fft_filter ! (ikr,ikz)
- double complex, DIMENSION(:,:), ALLOCATABLE :: fft_input ! (ikr,ikz)
-
- ! For WENO gradients
- real(dp), DIMENSION(:), ALLOCATABLE :: tmp1
- real(dp), DIMENSION(:), ALLOCATABLE :: tmp2
- real(dp), DIMENSION(:,:), ALLOCATABLE :: omegak_times_q3k
-
-
+ !COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE:: moments_Apl, moments_Bpl, moments_Cpl, moments_Dpl,&
+ ! moments_Epl, moments_Fpl, moments_Gpl, moments_Hpl
END MODULE array
diff --git a/src/auxval.F90 b/src/auxval.F90
index 5e6973ed18e89f1fcc64cd06b203f5564502b462..8925ac1c225e93843aa1a39e3517bed82f2882af 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -2,123 +2,18 @@ subroutine auxval
! Set auxiliary values, at beginning of simulation
USE basic
- USE space_grid
+ USE fourier_grid
USE array
- USE model, ONLY: gradient_scheme, diff_theta, diff_temp, diff_vpar, diff_moments
- USE gradients
! use mumps_bsplines, only: putrow_mumps => putrow, updtmat_mumps => updtmat, factor_mumps => factor, bsolve_mumps => bsolve ! from BSPLINES
- use fft
use prec_const
IMPLICIT NONE
INTEGER :: irows,irowe, irow, icol
- real(dp), DIMENSION(-2:2) :: stencil ! For laplace operator (needed for poisson solver)
-
WRITE(*,'(a/)') '=== Set auxiliary values ==='
- call set_pgrid
- call set_zgrid
-
-
- neq = (ize-izs+1) + 1 ! Number of equations to solve for poisson = Nz + 1 (additional equation to fix the constant)
+ call set_krgrid
+ call set_kzgrid
+ call set_pj
CALL memory ! Allocate memory for global arrays
-
- CALL set_gradient_scheme(gradient_scheme) ! Point to correct gradient functions
- if ((gradient_scheme .eq. 'up2') .and. ((diff_theta .ne. 0._dp) .or. (diff_temp .ne. 0._dp) .or. (diff_vpar .ne. 0._dp) .or. (diff_moments .ne. 0._dp))) then
- write(*,*) 'ERROR: the scheme is upwind but the numerical diffusion coefficients are non-zero. Program ends.'
- nlend = .true.
- endif
-
-
-! Setting up laplace matrix, eg (-2,1,0,0; 1,-2,1,0; 0,1,-2,1; 0,0,1,-2), to solve poisson equation
-
-! call set_laplace_stencil(gradient_scheme, stencil) ! for fa2 it is (0,1,-2,1,0)/dx**2
- ! call set_gradient_stencil(gradient_scheme, stencil) ! for fa2 it is (0,-1/2,0,1/2,0)/dx
-
-! irows = 1
-! irowe = neq-1
-
-! do irow=irows,irowe
-! do icol=irows,irowe
- ! do isten=max(-2,irows-irow),min(2,irowe-irow) ! loop over stencil
-! if ( (icol-irow .le. 2) .and. (icol-irow .ge. -2) ) then
-! call putele(laplace_smumps, irow, icol, stencil(icol-irow))
-! else
-! call putele(laplace_smumps, irow, icol, 0._dp)
-! endif
-! enddo
-! enddo
-
- ! Periodic boundary conditions. Eg transform last row from (0,0,0,1,-2) to (1,0,0,1,-2)
-! call putele(laplace_smumps, irows, irowe-1, stencil(-2))
-! call putele(laplace_smumps, irows, irowe, stencil(-1))
-
-! call putele(laplace_smumps, irows+1, irowe, stencil(-2))
-
-! call putele(laplace_smumps, irowe-1, irows, stencil(2))
-
-! call putele(laplace_smumps, irowe, irows, stencil(1))
-! call putele(laplace_smumps, irowe, irows+1, stencil(2))
-
-
- ! ! Boundary condition to fix the value of elec at position iz=21
- ! do icol=irows,irowe
- ! call putele(laplace_smumps, irows+20, icol, 0._dp)
- ! enddo
- ! call putele(laplace_smumps, irows+20, irows+20, 1._dp)
-
-
-
- ! ! Boundary condition to fix the mean value of electric field (replaces equation for iz=21)
- ! do icol=irows,irowe
- ! call putele(laplace_smumps, irows+20, icol, 1._dp/real(ize-izs+1,dp) )
- ! enddo
-
-
- ! Lagrange multiplier, adds an equation at neq^th row to fix the constant of integration by fixing the average value
-! do icol=irows,irowe
-! call putele(laplace_smumps, icol, neq, 1._dp)
-! call putele(laplace_smumps, neq, icol, deltaz)
-! enddo
-! call putele(laplace_smumps, neq, neq, 0._dp)
-
-
-
-
-! call factor(laplace_smumps) ! Factor
-
-
-
- CALL initialize_filter
-
-
-
-
-CONTAINS
-
- subroutine set_laplace_stencil(gradient_scheme, laplace_stencil)
-
- character(3), intent(in) :: gradient_scheme
- real(dp), dimension(-2:2), intent(out) :: laplace_stencil
-
- laplace_stencil = 0.0_dp
-
- select case (gradient_scheme)
- ! case('fa2','up2')
- ! laplace_stencil(-1) = 1._dp * deltazisq
- ! laplace_stencil(0) = -2._dp * deltazisq
- ! laplace_stencil(1) = 1._dp * deltazisq
-
- case('fd4','we4')
- laplace_stencil(-2) = -deltazisq/12._dp
- laplace_stencil(-1) = deltazisq*16._dp/12._dp
- laplace_stencil(0) = -deltazisq*30._dp/12._dp
- laplace_stencil(1) = deltazisq*16._dp/12._dp
- laplace_stencil(2) = -deltazisq/12._dp
-
- end select
-
- end subroutine set_laplace_stencil
-
END SUBROUTINE auxval
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index c4cc0ffc28a3050d049185e79935ce4b4577f057..0dcbedafd5018ffdc1d211685e7ea89a8af6f652 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -9,18 +9,15 @@ MODULE basic
real(dp) :: dt=1.0 ! Time step
real(dp) :: time=0 ! Current simulation time (Init from restart file)
-
INTEGER :: jobnum=0 ! Job number
INTEGER :: step=0 ! Calculation step of this run
INTEGER :: cstep=0 ! Current step number (Init from restart file)
LOGICAL :: nlend=.FALSE. ! Signal end of run
-
INTEGER :: ierr ! flag for MPI error
INTEGER :: iframe1d ! counting the number of times 1d datasets are outputed (for diagnose)
INTEGER :: iframe2d ! counting the number of times 2d datasets are outputed (for diagnose)
-
-
+ INTEGER :: iframe3d ! counting the number of times 3d datasets are outputed (for diagnose)
! List of logical file units
INTEGER :: lu_in = 90 ! File duplicated from STDIN
@@ -75,24 +72,28 @@ CONTAINS
ALLOCATE(a(is1:ie1))
a=0.0_dp
END SUBROUTINE allocate_array_dp1
+
SUBROUTINE allocate_array_dp2(a,is1,ie1,is2,ie2)
real(dp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
ALLOCATE(a(is1:ie1,is2:ie2))
a=0.0_dp
END SUBROUTINE allocate_array_dp2
+
SUBROUTINE allocate_array_dp3(a,is1,ie1,is2,ie2,is3,ie3)
real(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3))
a=0.0_dp
END SUBROUTINE allocate_array_dp3
+
SUBROUTINE allocate_array_dp4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
real(dp), DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4))
a=0.0_dp
END SUBROUTINE allocate_array_dp4
+
SUBROUTINE allocate_array_dp5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
@@ -100,30 +101,36 @@ CONTAINS
a=0.0_dp
END SUBROUTINE allocate_array_dp5
+ !========================================
+
SUBROUTINE allocate_array_dc1(a,is1,ie1)
DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1
ALLOCATE(a(is1:ie1))
a=CMPLX(0.0_dp,0.0_dp)
END SUBROUTINE allocate_array_dc1
+
SUBROUTINE allocate_array_dc2(a,is1,ie1,is2,ie2)
DOUBLE COMPLEX, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
ALLOCATE(a(is1:ie1,is2:ie2))
a=CMPLX(0.0_dp,0.0_dp)
END SUBROUTINE allocate_array_dc2
+
SUBROUTINE allocate_array_dc3(a,is1,ie1,is2,ie2,is3,ie3)
DOUBLE COMPLEX, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3))
a=CMPLX(0.0_dp,0.0_dp)
END SUBROUTINE allocate_array_dc3
+
SUBROUTINE allocate_array_dc4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
DOUBLE COMPLEX, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4))
a=CMPLX(0.0_dp,0.0_dp)
END SUBROUTINE allocate_array_dc4
+
SUBROUTINE allocate_array_dc5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
@@ -131,30 +138,36 @@ CONTAINS
a=CMPLX(0.0_dp,0.0_dp)
END SUBROUTINE allocate_array_dc5
+ !========================================
+
SUBROUTINE allocate_array_i1(a,is1,ie1)
INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1
ALLOCATE(a(is1:ie1))
a=0
END SUBROUTINE allocate_array_i1
+
SUBROUTINE allocate_array_i2(a,is1,ie1,is2,ie2)
INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
ALLOCATE(a(is1:ie1,is2:ie2))
a=0
END SUBROUTINE allocate_array_i2
+
SUBROUTINE allocate_array_i3(a,is1,ie1,is2,ie2,is3,ie3)
INTEGER, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3))
a=0
END SUBROUTINE allocate_array_i3
+
SUBROUTINE allocate_array_i4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4))
a=0
END SUBROUTINE allocate_array_i4
+
SUBROUTINE allocate_array_i5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
@@ -162,34 +175,41 @@ CONTAINS
a=0
END SUBROUTINE allocate_array_i5
+ !========================================
+
SUBROUTINE allocate_array_l1(a,is1,ie1)
LOGICAL, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1
ALLOCATE(a(is1:ie1))
a=.false.
END SUBROUTINE allocate_array_l1
+
SUBROUTINE allocate_array_l2(a,is1,ie1,is2,ie2)
LOGICAL, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
ALLOCATE(a(is1:ie1,is2:ie2))
a=.false.
END SUBROUTINE allocate_array_l2
+
SUBROUTINE allocate_array_l3(a,is1,ie1,is2,ie2,is3,ie3)
LOGICAL, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3))
a=.false.
END SUBROUTINE allocate_array_l3
+
SUBROUTINE allocate_array_l4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4))
a=.false.
END SUBROUTINE allocate_array_l4
+
SUBROUTINE allocate_array_l5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5))
a=.false.
END SUBROUTINE allocate_array_l5
+
END MODULE basic
diff --git a/src/control.F90 b/src/control.F90
index ebf206747dd04761f80e6339733d20cb53283bea..e2d481ec18e950f6b9a2d179b70f48c8a846f54b 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -34,7 +34,7 @@ SUBROUTINE control
!
! 1.5 Initial conditions
WRITE(*,*) 'Create initial state...'
- CALL inital
+ CALL inital
WRITE(*,*) '...initial state created'
! 1.6 Initial diagnostics
@@ -46,18 +46,18 @@ SUBROUTINE control
!________________________________________________________________________________
! 2. Main loop
- DO
- step = step+1
- cstep = cstep+1
- CALL stepon
- time = time + dt
+ !DO
+ ! step = step+1
+ ! cstep = cstep+1
+ ! CALL stepon
+ ! time = time + dt
- CALL tesend
- CALL diagnose(step)
- IF( nlend ) EXIT ! exit do loop
+! CALL tesend
+! CALL diagnose(step)
+! IF( nlend ) EXIT ! exit do loop
! CALL write_restart ! if want to write a restart file every so often (in case of crash)
- END DO
+! END DO
!________________________________________________________________________________
! 9. Epilogue
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 6664db10327003b157501299a61b17a249f80c9e..ae0ef8575e4b661a5ee4c4f11213a2f2ab75ca09 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -2,7 +2,7 @@ SUBROUTINE diagnose(kstep)
! Diagnostics, writing simulation state to disk
USE basic
- USE space_grid
+ USE fourier_grid
USE diagnostics_par
USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach
USE model
@@ -47,68 +47,56 @@ SUBROUTINE diagnose(kstep)
! Data group
CALL creatg(fidres, "/data", "data")
- ! CALL creatg(fidres, "/data/var0d", "0d history arrays")
- CALL creatg(fidres, "/data/var1d", "1d profiles")
+ !CALL creatg(fidres, "/data/var0d", "0d history arrays")
+ !CALL creatg(fidres, "/data/var1d", "1d profiles")
CALL creatg(fidres, "/data/var2d", "2d profiles")
- ! CALL creatg(fidres, "/data/var3d", "3d profiles")
+ CALL creatg(fidres, "/data/var3d", "3d profiles")
! Initialize counter of number of saves for each category
- IF (cstep==0) THEN
- iframe1d=0
- END IF
- CALL attach(fidres,"/data/var1d/" , "frames", iframe1d)
+ !IF (cstep==0) THEN
+ ! iframe1d=0
+ !END IF
+ !CALL attach(fidres,"/data/var1d/" , "frames", iframe1d)
IF (cstep==0) THEN
iframe2d=0
END IF
CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
+ IF (cstep==0) THEN
+ iframe3d=0
+ END IF
+ CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
! File group
CALL creatg(fidres, "/files", "files")
CALL attach(fidres, "/files", "jobnum", jobnum)
-
- ! var0d group (empty in our case)
+ ! var2d group
rank = 0
-! dims(1) = 0
- ! var1d group
- CALL creatd(fidres, rank, dims, "/data/var1d/time", "Time t*w_pe") ! time of the saves for 1d variables
- CALL creatd(fidres, rank, dims, "/data/var1d/cstep", "iteration number") ! cstep of the saves for 1d variables
-
- IF (write_theta) THEN
- CALL creatg(fidres, "/data/var1d/theta", "theta")
- CALL putarr(fidres, "/data/var1d/theta/coordz", zarray(izs:ize), "z/lambda_e0",ionode=0)
- END IF
- IF (write_temp) THEN
- CALL creatg(fidres, "/data/var1d/temp", "temp")
- CALL putarr(fidres, "/data/var1d/temp/coordz", zarray(izs:ize), "z/lambda_e0",ionode=0)
- END IF
- IF (write_vpar) THEN
- CALL creatg(fidres, "/data/var1d/vpar", "vpar")
- CALL putarr(fidres, "/data/var1d/vpar/coordz", zarray(izs:ize), "z/lambda_e0",ionode=0)
- END IF
+ CALL creatd(fidres, rank, dims, "/data/var2d/time", "Time t*c_s/R")
+ CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number")
IF (write_phi) THEN
- CALL creatg(fidres, "/data/var1d/phi", "phi")
- CALL putarr(fidres, "/data/var1d/phi/coordz", zarray(izs:ize), "z/lambda_e0",ionode=0)
+ CALL creatg(fidres, "/data/var2d/phi", "phi")
+ CALL putarr(fidres, "/data/var2d/phi/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var2d/phi/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
END IF
-
- ! var2d group
+ ! var3d group
rank = 0
- CALL creatd(fidres, rank, dims, "/data/var2d/time", "Time t*c_s/R")
- CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number")
-
+ CALL creatd(fidres, rank, dims, "/data/var3d/time", "Time t*c_s/R")
+ CALL creatd(fidres, rank, dims, "/data/var3d/cstep", "iteration number")
IF (write_moments) THEN
- CALL creatg(fidres, "/data/var2d/moments", "moments")
- CALL putarr(fidres, "/data/var2d/moments/coordp", parray(ips:ipe), "p",ionode=0)
- CALL putarr(fidres, "/data/var2d/moments/coordz", zarray(izs:ize), "z/lambda_e0",ionode=0)
+ CALL creatg(fidres, "/data/var3d/moments", "moments")
+ CALL putarr(fidres, "/data/var3d/moments/coordpj", pjarray(ipjs:ipje),"(Jmaxa+1)*p+j+1",ionode=0)
+ CALL putarr(fidres, "/data/var3d/moments/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var3d/moments/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
END IF
! Add input namelist variables as attributes of /data/input, defined in srcinfo.h
- WRITE(*,*) 'VERSION=',VERSION
- WRITE(*,*) 'BRANCH=',BRANCH
- WRITE(*,*) 'AUTHOR=',AUTHOR
- WRITE(*,*) 'HOST=',HOST
+ WRITE(*,*) 'VERSION=', VERSION
+ WRITE(*,*) 'BRANCH=', BRANCH
+ WRITE(*,*) 'AUTHOR=', AUTHOR
+ WRITE(*,*) 'HOST=', HOST
IF(jobnum .LE. 99) THEN
WRITE(str,'(a,i2.2)') "/data/input.",jobnum
@@ -117,18 +105,18 @@ SUBROUTINE diagnose(kstep)
END IF
rank=0
CALL creatd(fidres, rank,dims,TRIM(str),'Input parameters')
- CALL attach(fidres, TRIM(str), "version", VERSION) !defined in srcinfo.h
- CALL attach(fidres, TRIM(str), "branch", BRANCH) !defined in srcinfo.h
- CALL attach(fidres, TRIM(str), "author", AUTHOR) !defined in srcinfo.h
- CALL attach(fidres, TRIM(str), "execdate", EXECDATE) !defined in srcinfo.h
- CALL attach(fidres, TRIM(str), "host", HOST) !defined in srcinfo.h
- CALL attach(fidres, TRIM(str), "start_time", time)
- CALL attach(fidres, TRIM(str), "start_cstep", cstep)
- CALL attach(fidres, TRIM(str), "dt", dt)
- CALL attach(fidres, TRIM(str), "tmax", tmax)
- CALL attach(fidres, TRIM(str), "nrun", nrun)
-
- CALL space_grid_outputinputs(fidres, str)
+ CALL attach(fidres, TRIM(str), "version", VERSION) !defined in srcinfo.h
+ CALL attach(fidres, TRIM(str), "branch", BRANCH) !defined in srcinfo.h
+ CALL attach(fidres, TRIM(str), "author", AUTHOR) !defined in srcinfo.h
+ CALL attach(fidres, TRIM(str), "execdate", EXECDATE) !defined in srcinfo.h
+ CALL attach(fidres, TRIM(str), "host", HOST) !defined in srcinfo.h
+ CALL attach(fidres, TRIM(str), "start_time", time)
+ CALL attach(fidres, TRIM(str), "start_cstep", cstep)
+ CALL attach(fidres, TRIM(str), "dt", dt)
+ CALL attach(fidres, TRIM(str), "tmax", tmax)
+ CALL attach(fidres, TRIM(str), "nrun", nrun)
+
+ CALL fourier_grid_outputinputs(fidres, str)
CALL output_par_outputinputs(fidres, str)
@@ -167,11 +155,7 @@ SUBROUTINE diagnose(kstep)
END IF
! 2.2 1d profiles
- IF (nsave_1d .NE. 0) THEN
- IF ( MOD(cstep, nsave_1d) == 0 ) THEN
- CALL diagnose_1d
- END IF
- END IF
+ ! empty in our case
! 2.3 2d profiles
IF (nsave_2d .NE. 0) THEN
@@ -181,8 +165,11 @@ SUBROUTINE diagnose(kstep)
END IF
! 2.4 3d profiles
- ! empty in our case
-
+ IF (nsave_3d .NE. 0) THEN
+ IF (MOD(cstep, nsave_3d) == 0) THEN
+ CALL diagnose_3d
+ END IF
+ END IF
!________________________________________________________________________________
! 3. Final diagnostics
@@ -196,9 +183,6 @@ SUBROUTINE diagnose(kstep)
END SUBROUTINE diagnose
-
-
-
SUBROUTINE diagnose_0d
USE basic
@@ -218,63 +202,6 @@ SUBROUTINE diagnose_0d
!CALL append(fidres,"/data/var0d/cstep" ,real(cstep,dp), ionode=0)
END SUBROUTINE diagnose_0d
-
-
-
-SUBROUTINE diagnose_1d
-
- USE basic
- USE futils, ONLY: append, getatt, attach, putarr
- USE fields
- USE time_integration
- USE diagnostics_par
- use prec_const
- IMPLICIT NONE
-
-
- CALL append(fidres, "/data/var1d/time", time ,ionode=0)
- CALL append(fidres, "/data/var1d/cstep", real(cstep,dp),ionode=0)
- CALL getatt(fidres,"/data/var1d/" , "frames", iframe1d)
- iframe1d=iframe1d+1
- CALL attach(fidres,"/data/var1d/" , "frames", iframe1d)
-
- IF (write_theta) THEN
- CALL write_field1d(theta(:, updatetlevel), 'theta')
- END IF
- IF (write_temp) THEN
- CALL write_field1d(temp(:, updatetlevel), 'temp')
- END IF
- IF (write_vpar) THEN
- CALL write_field1d(vpar(:, updatetlevel), 'vpar')
- END IF
- IF (write_phi) THEN
- CALL write_field1d(phi(:), 'phi')
- END IF
-CONTAINS
-
- SUBROUTINE write_field1d(field, text)
- USE futils, ONLY: attach, putarr
- USE space_grid, only: izs, ize
- use prec_const
- IMPLICIT NONE
-
- real(dp), DIMENSION(izs:ize), INTENT(IN) :: field
- CHARACTER(*), INTENT(IN) :: text
-
- CHARACTER(LEN=50) :: dset_name
-
-
- WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var1d", TRIM(text), iframe1d
- CALL putarr(fidres, dset_name, field(izs:ize), text, ionode=0)
-
-
- CALL attach(fidres, dset_name, "time", time)
-
- END SUBROUTINE write_field1d
-
-END SUBROUTINE diagnose_1d
-
-
SUBROUTINE diagnose_2d
@@ -287,31 +214,31 @@ SUBROUTINE diagnose_2d
use prec_const
IMPLICIT NONE
- CALL append(fidres, "/data/var2d/time", time ,ionode=0)
- CALL append(fidres, "/data/var2d/cstep", real(cstep,dp),ionode=0)
- CALL getatt(fidres,"/data/var2d/" , "frames", iframe2d)
+ CALL append(fidres, "/data/var2d/time", time,ionode=0)
+ CALL append(fidres, "/data/var2d/cstep", real(cstep,dp),ionode=0)
+ CALL getatt(fidres, "/data/var2d/", "frames",iframe2d)
iframe2d=iframe2d+1
CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
- IF (write_moments) THEN
- CALL write_field2d(moments(:,:,updatetlevel), 'moments')
+ IF (write_phi) THEN
+ CALL write_field2d(phi(:,:), 'phi')
END IF
CONTAINS
SUBROUTINE write_field2d(field, text)
USE futils, ONLY: attach, putarr
- USE space_grid, only: ips, ipe, izs, ize
+ USE fourier_grid, only: ikrs,ikre, ikzs,ikze
use prec_const
IMPLICIT NONE
- real(dp), DIMENSION(ips:ipe,izs:ize), INTENT(IN) :: field
+ COMPLEX(dp), DIMENSION(ikrs:ikre, ikzs:ikze), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
CHARACTER(LEN=50) :: dset_name
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var2d", TRIM(text), iframe2d
- CALL putarr(fidres, dset_name, field(ips:ipe,izs:ize),ionode=0)
+ CALL putarr(fidres, dset_name, field(ikrs:ikre, ikzs:ikze),ionode=0)
CALL attach(fidres, dset_name, "time", time)
@@ -319,3 +246,44 @@ CONTAINS
END SUBROUTINE diagnose_2d
+SUBROUTINE diagnose_3d
+
+ USE basic
+ USE futils, ONLY: append, getatt, attach, putarrnd
+ USE fields
+ USE time_integration
+ USE diagnostics_par
+ use prec_const
+ IMPLICIT NONE
+
+ CALL append(fidres, "/data/var3d/time", time,ionode=0)
+ CALL append(fidres, "/data/var3d/cstep", real(cstep,dp),ionode=0)
+ CALL getatt(fidres, "/data/var3d/", "frames",iframe3d)
+ iframe3d=iframe3d+1
+ CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
+
+ IF (write_moments) THEN
+ CALL write_field3d(moments(:,:,:,updatetlevel), 'moments')
+ END IF
+
+ CONTAINS
+
+ SUBROUTINE write_field3d(field, text)
+ USE futils, ONLY: attach, putarr
+ USE fourier_grid, only: ipjs,ipje, ikrs,ikre, ikzs,ikze
+ use prec_const
+ IMPLICIT NONE
+
+ COMPLEX(dp), DIMENSION(ipjs:ipje,ikrs:ikre,ikzs:ikze), INTENT(IN) :: field
+ CHARACTER(*), INTENT(IN) :: text
+
+ CHARACTER(LEN=50) :: dset_name
+
+ WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
+ CALL putarr(fidres, dset_name, field(ipjs:ipje,ikrs:ikre,ikzs:ikze),ionode=0)
+
+ CALL attach(fidres, dset_name, "time", time)
+
+ END SUBROUTINE write_field3d
+
+ END SUBROUTINE diagnose_3d
\ No newline at end of file
diff --git a/src/fft_mod.F90 b/src/fft_mod.F90
deleted file mode 100644
index 2aa57198c01e1832b6885b507c79a235a65f1689..0000000000000000000000000000000000000000
--- a/src/fft_mod.F90
+++ /dev/null
@@ -1,99 +0,0 @@
-module fft
- ! Fast fourier transform module
- ! Used to supress high frequency instabilities
-
- use prec_const
- use array
- implicit none
- private
-
-
- ! Public Functions
- PUBLIC :: fft_recursive, initialize_filter, suppress_high_freq
-
-
- contains
-
- ! In place Cooley-Tukey FFT, tsource : rosettacode.org/wiki/Fast_Fourier_transform
- recursive subroutine fft_recursive(x)
- complex(kind=dp), dimension(:), intent(inout) :: x
- complex(kind=dp) :: t
- integer :: N
- integer :: i
- complex(kind=dp), dimension(:), allocatable :: even, odd
-
- N=size(x)
-
- if(N .le. 1) return
-
- allocate(odd((N+1)/2))
- allocate(even(N/2))
-
- ! divide
- odd =x(1:N:2)
- even=x(2:N:2)
-
- ! conquer
- call fft_recursive(odd)
- call fft_recursive(even)
-
- ! combine
- do i=1,N/2
- t=exp(cmplx(0.0_dp,-2.0_dp*pi*real(i-1,dp)/real(N,dp),kind=dp))*even(i)
- x(i) = odd(i) + t
- x(i+N/2) = odd(i) - t
- end do
-
- deallocate(odd)
- deallocate(even)
- end subroutine fft_recursive
-
-
- subroutine initialize_filter
- use space_grid
- use array
- use model
- implicit none
- integer :: i
- integer :: N
- integer :: center
-
- N = (ize-izs+1)
- center = int( ceiling( (N+1)*0.5 ) )
-
- do i=izs,ize
- fft_filter(i) = 1._dp- exp( -1._dp*(i-center)*(i-center)/ (2*fft_sigma*fft_sigma) ) ! or just a step function?
- ! write(*,'("(", F20.15, ",", F20.15, "i )")') fft_filter(i)
- end do
- end subroutine initialize_filter
-
-
- subroutine suppress_high_freq(x)
- use space_grid
- use array
- implicit none
- integer :: i
- real(dp), dimension(izs:ize), intent(inout) :: x
- real(dp) :: invN
-
- do i=izs,ize
- fft_input(i) = x(i)
- end do
-
- call fft_recursive(fft_input) ! Fast Fourier Transform
-
- do i=izs,ize ! Apply filter
- fft_input(i) = fft_input(i)*fft_filter(i)
- end do
-
- call fft_recursive(fft_input) ! Invert the FFT
-
- invN = 1._dp/(ize-izs+1._dp)
- do i=izs,ize
- x(i) = real(fft_input(ize-i+1))*invN
- end do
- end subroutine suppress_high_freq
-
-
-end module fft
-
diff --git a/src/fields_mod.F90 b/src/fields_mod.F90
index 91f75f765d740a7a78fb72339339cb46faafc192..00f8271986bfd0e4e88e975203d98144bed79e33 100644
--- a/src/fields_mod.F90
+++ b/src/fields_mod.F90
@@ -2,66 +2,14 @@ MODULE fields
use prec_const
implicit none
-
-!------------------DENSITY Na00------------------
-
- ! Natural logarithm of normalized density: \ln \hat{N}
-! real(dp), DIMENSION(:,:), ALLOCATABLE :: theta
-
- ! Natural logarithm of normalized electrons density: \ln \hat{N}
- real(dp), DIMENSION(:,:,:), ALLOCATABLE :: theta_i ! (kr,kz,updatetlevel)
-
- ! Natural logarithm of normalized ions density: \ln \hat{N}
- real(dp), DIMENSION(:,:,:), ALLOCATABLE :: theta_e ! (kr,kz,updatetlevel)
-
-!------------------FLUID VELOCITY Na10------------------
-
- ! Normalized electron gyrocenter parallel fluid velocity \hat{u}
-! real(dp), DIMENSION(:,:), ALLOCATABLE :: vpar
-
- ! Normalized electron gyrocenter parallel fluid velocity \hat{u}
- real(dp), DIMENSION(:,:,:), ALLOCATABLE :: v_par_e ! (kr,kz,updatetlevel)
-
- ! Normalized ion gyrocenter parallel fluid velocity \hat{u}
- real(dp), DIMENSION(:,:,:), ALLOCATABLE :: v_par_i ! (kr,kz,updatetlevel)
-
-!------------------PARALLEL TEMPERATURE Na20------------------
-
- ! Natural logarithm of normalized parallel electron temperature: \ln \hat{T}
-! real(dp), DIMENSION(:,:,:), ALLOCATABLE :: temp
-
- ! Natural logarithm of normalized parallel electron temperature: \ln \hat{T}
- real(dp), DIMENSION(:,:,:), ALLOCATABLE :: temp_par_e ! (kr,kz,updatetlevel)
-
- ! Natural logarithm of normalized perpendicular electron temperature: \ln \hat{T}
- real(dp), DIMENSION(:,:,:), ALLOCATABLE :: temp_par_i ! (kr,kz,updatetlevel)
-
-!------------------PERPENDICULAR TEMPERATURE Na01------------------
-
- ! Natural logarithm of normalized parallel ion temperature: \ln \hat{T}
- real(dp), DIMENSION(:,:,:), ALLOCATABLE :: temp_per_e ! (kr,kz,updatetlevel)
-
- ! Natural logarithm of normalized perpendicular ion temperature: \ln \hat{T}
- real(dp), DIMENSION(:,:,:), ALLOCATABLE :: temp_per_i ! (kr,kz,updatetlevel)
-
!------------------MOMENTS Napj------------------
-
- ! Hermite-Moments: N_e^pj ! dimensions correspond to: moments p (from 3 to pmax), z, updatetlevel.
-! real(dp), DIMENSION(:,:,:), ALLOCATABLE :: moments
-
- ! Hermite-Moments: N_e^pj ! dimensions correspond to: moments pj (from 3 to pmax, 1 to jmax), kr, kz, updatetlevel.
- real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_e
-
- ! Hermite-Moments: N_i^pj ! dimensions correspond to: moments pj (from 3 to pmax, 1 to jmax), kr, kz, updatetlevel.
- real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_i
+ ! Hermite-Moments: N_a^pj ! dimensions correspond to: ipj, kr, kz, updatetlevel.
+ COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments
!------------------ELECTROSTATIC POTENTIAL------------------
- ! Normalized electric potential: \hat{\phi}
-! real(dp), DIMENSION(:,:), ALLOCATABLE :: phi
-
- ! Normalized electric potential: \hat{\phi}
- real(dp), DIMENSION(:,:,:), ALLOCATABLE :: phi ! (kr,kz,updatetlevel)
+ ! Normalized electric potential: \hat{\phi} ! (kr,kz)
+ COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: phi
END MODULE fields
diff --git a/src/fourier_grid_mod.F90 b/src/fourier_grid_mod.F90
index 94ea80fde7d496b35cdd855ae02d172002558610..b4d7232070bc12058d8e042ac61d3096767a4f67 100644
--- a/src/fourier_grid_mod.F90
+++ b/src/fourier_grid_mod.F90
@@ -6,86 +6,91 @@ MODULE fourier_grid
PRIVATE
! GRID Namelist
- INTEGER, PUBLIC, PROTECTED :: pmax=3 ! The maximal Hermite-moment computed
- INTEGER, PUBLIC, PROTECTED :: jmax=2 ! The maximal Laguerre-moment computed
- INTEGER, PUBLIC, PROTECTED :: nkr=1 ! Number of total internal grid points in kr
- INTEGER, PUBLIC, PROTECTED :: nkz=1 ! Number of total internal grid points in kz
- REAL(dp), PUBLIC, PROTECTED :: krmin=0._dp ! kr coordinate for left boundary
- REAL(dp), PUBLIC, PROTECTED :: krmax=1._dp ! kr coordinate for right boundary
- REAL(dp), PUBLIC, PROTECTED :: kzmin=0._dp ! kz coordinate for left boundary
- REAL(dp), PUBLIC, PROTECTED :: kzmax=1._dp ! kz coordinate for right boundary
-
+ INTEGER, PUBLIC, PROTECTED :: pmaxe = 2 ! The maximal electron Hermite-moment computed
+ INTEGER, PUBLIC, PROTECTED :: jmaxe = 2 ! The maximal electron Laguerre-moment computed
+ INTEGER, PUBLIC, PROTECTED :: pmaxi = 2 ! The maximal ion Hermite-moment computed
+ INTEGER, PUBLIC, PROTECTED :: jmaxi = 2 ! The maximal ion Laguerre-moment computed
+ INTEGER, PUBLIC, PROTECTED :: nkr = 10 ! Number of total internal grid points in kr
+ REAL(dp), PUBLIC, PROTECTED :: krmin = 0._dp ! kr coordinate for left boundary
+ REAL(dp), PUBLIC, PROTECTED :: krmax = 1._dp ! kr coordinate for right boundary
+ INTEGER, PUBLIC, PROTECTED :: nkz = 10 ! Number of total internal grid points in kz
+ REAL(dp), PUBLIC, PROTECTED :: kzmin = 0._dp ! kz coordinate for left boundary
+ REAL(dp), PUBLIC, PROTECTED :: kzmax = 1._dp ! kz coordinate for right boundary
! Indices of s -> start , e-> end
- INTEGER, PUBLIC, PROTECTED :: ips, ipe, ijs, ije
- INTEGER, PUBLIC, PROTECTED :: ikrs, ikre, ikzs, ikze
+ INTEGER, PUBLIC, PROTECTED :: ipjs, ipje
+ INTEGER, PUBLIC, PROTECTED :: Nmome, Nmomi, Nmomtot
+ INTEGER, PUBLIC, PROTECTED :: ikrs, ikre, ikzs, ikze
! Toroidal direction
- REAL(dp), PUBLIC, PROTECTED :: deltakr, deltakz
- real(dp), PUBLIC, PROTECTED :: deltakri, deltakzi
- real(dp), PUBLIC, PROTECTED :: deltakrih, deltakzih
- real(dp), PUBLIC, PROTECTED :: deltakrisq, deltakzisq
-
- ! Grids containing position
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: parray
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: zarray
-
- INTEGER, PUBLIC :: neq ! Number of equations in poisson solver (number of rows/cols of laplace matrix)
+ REAL(dp), PUBLIC, PROTECTED :: deltakr
+ REAL(dp), PUBLIC, PROTECTED :: deltakz
+
+ ! Grids containing position in fourier space
+ REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: kzarray
+ REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: krarray
+ ! Grid containing the polynomials degrees
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: pjarray
+
! Public Functions
- PUBLIC :: set_pgrid, set_zgrid
- PUBLIC :: space_grid_readinputs, space_grid_outputinputs
+ PUBLIC :: set_krgrid, set_kzgrid, set_pj
+ PUBLIC :: fourier_grid_readinputs, fourier_grid_outputinputs
+ PUBLIC :: bare, bari
+ PUBLIC :: rabe, rabi
contains
- subroutine set_pjgrid
- ! Initialize pj grid (Hermite moment) parameters
-
+ subroutine set_krgrid
use prec_const
implicit none
- integer :: ip, ij
- !since N_e^00=1 and N_e^10 = N_e^20 = = N_e^01 = 0
- ips=0
- ijs=0
- ipe=pmax
- jpe=jmax
-
- allocate(pjarray(ips:ipe,jps:jpe))
- DO ip = ips,ipe
- DO ij = ijs,ije
- pjarray(ip,ij) = ip ! simply contains the moment number
- END DO
+ integer :: ikr
+ ! Start and end indices of grid
+ ikrs = 1
+ ikre = nkr
+ ! Grid spacings, precompute some inverses
+ deltakr = (krmax - krmin) / real(nkr,dp)
+ ! Discretized kr positions
+ allocate(krarray(ikrs:ikre))
+ DO ikr = ikrs,ikre
+ krarray(ikr) = krmin + real(ikr-1,dp) * deltakr
END DO
+ end subroutine set_krgrid
- end subroutine set_pjgrid
-
-
- subroutine set_zgrid
- ! Initialize z grid
-
+ subroutine set_kzgrid
use prec_const
implicit none
-
- integer :: iz
-
+ integer :: ikz
! Start and end indices of grid
- izs=1
- ize=nz
-
+ ikzs = 1
+ ikze = nkz
! Grid spacings, precompute some inverses
- deltaz = (zmax-zmin)/real(nz,dp) !zmax*2._dp*pi/real(nz,dp)
- deltazi = 1._dp/deltaz ! inverse
- deltazih = 0.5_dp/deltaz ! half of inverse
- deltazisq = 1._dp/deltaz/deltaz ! inverse squared
-
- ! Discretized z positions
- allocate(zarray(izs:ize))
- DO iz = izs,ize
- zarray(iz) = zmin + real(iz-1,dp)*deltaz
+ deltakz = (kzmax - kzmin) / real(nkz,dp)
+ ! Discretized kz positions
+ allocate(kzarray(ikzs:ikze))
+ DO ikz = ikzs,ikze
+ kzarray(ikz) = kzmin + real(ikz-1,dp) * deltakz
END DO
+ end subroutine set_kzgrid
- end subroutine set_zgrid
-
+ subroutine set_pj
+ use prec_const
+ implicit none
+ integer :: ipj
+ ! number of electrons moments
+ Nmome = (Pmaxe + 1)*(Jmaxe + 1)
+ ! number of ions moments
+ Nmomi = (Pmaxi + 1)*(Jmaxi + 1)
+ ! total number of moments
+ Nmomtot = Nmome + Nmomi
+ ipjs = 1
+ ipje = Nmomtot
+ ! Polynomials degrees pj = (Jmaxs + 1)*p + j + 1
+ allocate(pjarray(ipjs:ipje))
+ DO ipj = ipjs,ipje
+ pjarray(ipj) = ipj
+ END DO
+ end subroutine set_pj
SUBROUTINE fourier_grid_readinputs
! Read the input parameters
@@ -95,7 +100,8 @@ contains
USE prec_const
IMPLICIT NONE
- NAMELIST /GRID/ pmax, nz, zmin, zmax
+ NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
+ nkr, krmin, krmax, nkz, kzmin, kzmax
READ(lu_in,grid)
WRITE(*,grid)
@@ -113,13 +119,48 @@ contains
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
+ CALL attach(fidres, TRIM(str), "pmaxe", pmaxe)
+ CALL attach(fidres, TRIM(str), "jmaxe", jmaxe)
+ CALL attach(fidres, TRIM(str), "pmaxi", pmaxi)
+ CALL attach(fidres, TRIM(str), "jmaxi", jmaxi)
+ CALL attach(fidres, TRIM(str), "nkr", nkr)
+ CALL attach(fidres, TRIM(str), "krmin", krmin)
+ CALL attach(fidres, TRIM(str), "krmax", krmax)
+ CALL attach(fidres, TRIM(str), "nkz", nkz)
+ CALL attach(fidres, TRIM(str), "kzmin", kzmin)
+ CALL attach(fidres, TRIM(str), "kzmax", kzmax)
+ END SUBROUTINE fourier_grid_outputinputs
+
+ SUBROUTINE bare(p,j,idx)
+ USE prec_const
+ IMPLICIT NONE
+ INTEGER, INTENT(in) :: p,j
+ INTEGER, INTENT(out):: idx
+
+ idx = (jmaxe + 1)*p + j + 1
+
+ END SUBROUTINE bare
+
+ SUBROUTINE bari(p,j,idx)
+ INTEGER, INTENT(in) :: p,j
+ INTEGER, INTENT(out):: idx
+
+ idx = (jmaxi + 1)*p + j + 1
- CALL attach(fidres, TRIM(str), "pmax", pmax)
- CALL attach(fidres, TRIM(str), "nz", nz)
- CALL attach(fidres, TRIM(str), "zmin", zmin)
- CALL attach(fidres, TRIM(str), "zmax", zmax)
+ END SUBROUTINE bari
- END SUBROUTINE fourier_grid_outputinputs
+ SUBROUTINE rabe(idx, p, j)
+ INTEGER, INTENT(in) :: idx
+ INTEGER, INTENT(out):: p,j
+ p = FLOOR(real((idx-1) / (jmaxe + 1)))
+ j = idx - p * (jmaxe+1)
+ END SUBROUTINE rabe
+ SUBROUTINE rabi(idx, p, j)
+ INTEGER, INTENT(in):: idx
+ INTEGER, INTENT(out) :: p,j
+ p = FLOOR(real((idx-Nmome - 1) / (jmaxi + 1)))
+ j = (idx-Nmome) - p * (jmaxi+1)
+ END SUBROUTINE rabi
END MODULE fourier_grid
diff --git a/src/grad_fa2.h b/src/grad_fa2.h
deleted file mode 100644
index 799aeeefc8e62b8f6804c1ca0cddf316b2285862..0000000000000000000000000000000000000000
--- a/src/grad_fa2.h
+++ /dev/null
@@ -1,32 +0,0 @@
-SUBROUTINE gradz_fa2( f , f_z )
-
- use prec_const
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f
- real(dp), dimension(izs:ize),intent(out) :: f_z
- integer :: iz
-
- f_z(izs) = deltazih * ( f(izs+1) - f(ize) )
- do iz=izs+1,ize-1
- f_z(iz) = deltazih * ( f(iz+1) - f(iz-1) )
- end do
- f_z(ize) = deltazih * ( f(izs) - f(ize-1) )
-
-END SUBROUTINE gradz_fa2
-
-
-SUBROUTINE gradzz_fa2( f , f_z )
-
- use prec_const
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f
- real(dp), dimension(izs:ize),intent(out) :: f_z
- integer :: iz
-
- f_z(izs) = deltazisq * ( f(izs+1) - 2._dp*f(izs) + f(ize) )
- do iz=izs+1,ize-1
- f_z(iz) = deltazisq * ( f(iz+1) - 2._dp*f(iz) + f(iz-1) )
- end do
- f_z(ize) = deltazisq * ( f(izs) - 2._dp*f(ize) + f(ize-1) )
-
-END SUBROUTINE gradzz_fa2
diff --git a/src/grad_fd4.h b/src/grad_fd4.h
deleted file mode 100644
index 96b85887d1c96bf093a48749e2c150f2448dd378..0000000000000000000000000000000000000000
--- a/src/grad_fd4.h
+++ /dev/null
@@ -1,159 +0,0 @@
-SUBROUTINE gradz_fd4( f , f_z ) ! implemented for periodic boundary conditions
-
- use prec_const
- use space_grid
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f
- real(dp), dimension(izs:ize),intent(out) :: f_z
- real(dp), dimension(-2:2) :: coef ! fd4 stencil
- integer :: iz
-
- coef = (/ deltazi*onetwelfth, &
- -deltazi*2._dp*onethird, &
- 0._dp, &
- deltazi*2._dp*onethird, &
- -deltazi*onetwelfth /)
-
- f_z(izs) = coef(-2)*f(ize-1) + coef(-1)*f(ize) + coef(1)*f(izs+1) + coef(2)*f(izs+2)
- f_z(izs+1) = coef(-2)*f(ize) + coef(-1)*f(izs) + coef(1)*f(izs+2) + coef(2)*f(izs+3)
- do iz=izs+2,ize-2
- f_z(iz) = coef(-2)*f(iz-2) + coef(-1)*f(iz-1) + coef(1)*f(iz+1) + coef(2)*f(iz+2)
- end do
- f_z(ize-1) = coef(-2)*f(ize-3) + coef(-1)*f(ize-2) + coef(1)*f(ize) + coef(2)*f(izs)
- f_z(ize) = coef(-2)*f(ize-2) + coef(-1)*f(ize-1) + coef(1)*f(izs) + coef(2)*f(izs+1)
-
-END SUBROUTINE gradz_fd4
-
-
-SUBROUTINE gradz_n2v_fd4( f , f_z ) ! implemented for periodic boundary conditions
-
- use prec_const
- use space_grid
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f
- real(dp), dimension(izs:ize),intent(out) :: f_z
- real(dp), dimension(-2:1) :: coef ! fd4 stencil
- integer :: iz
-
- coef = (/ deltazi*onetwentyfourth, &
- -deltazi*nineeighths, &
- deltazi*nineeighths, &
- -deltazi*onetwentyfourth /)
-
- f_z(izs) = coef(-2)*f(ize-1) + coef(-1)*f(ize) + coef(0)*f(izs) + coef(1)*f(izs+1)
- f_z(izs+1) = coef(-2)*f(ize) + coef(-1)*f(izs) + coef(0)*f(izs+1) + coef(1)*f(izs+2)
- do iz=izs+2,ize-1
- f_z(iz) = coef(-2)*f(iz-2) + coef(-1)*f(iz-1) + coef(0)*f(iz) + coef(1)*f(iz+1)
- end do
- f_z(ize) = coef(-2)*f(ize-2) + coef(-1)*f(ize-1) + coef(0)*f(ize) + coef(1)*f(izs)
-
-END SUBROUTINE gradz_n2v_fd4
-
-
-SUBROUTINE gradz_v2n_fd4( f , f_z ) ! implemented for periodic boundary conditions
-
- use prec_const
- use space_grid
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f
- real(dp), dimension(izs:ize),intent(out) :: f_z
- real(dp), dimension(-1:2) :: coef ! fd4 stencil
- integer :: iz
-
- coef = (/ deltazi*onetwentyfourth, &
- -deltazi*nineeighths, &
- deltazi*nineeighths, &
- -deltazi*onetwentyfourth /)
-
- f_z(izs) = coef(-1)*f(ize) + coef(0)*f(izs) + coef(1)*f(izs+1) + coef(2)*f(izs+2)
- do iz=izs+1,ize-2
- f_z(iz) = coef(-1)*f(iz-1) + coef(0)*f(iz) + coef(1)*f(iz+1) + coef(2)*f(iz+2)
- end do
- f_z(ize-1) = coef(-1)*f(ize-2) + coef(0)*f(ize-1) + coef(1)*f(ize) + coef(2)*f(izs)
- f_z(ize) = coef(-1)*f(ize-1) + coef(0)*f(ize) + coef(1)*f(izs) + coef(2)*f(izs+1)
-
-END SUBROUTINE gradz_v2n_fd4
-
-
-
-SUBROUTINE gradzz_fd4( f , f_z ) ! implemented for periodic boundary conditions
-
- use prec_const
- use space_grid
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f
- real(dp), dimension(izs:ize),intent(out) :: f_z
- real(dp), dimension(-2:2) :: coef ! fd4 stencil
- integer :: iz
-
- coef = (/ -deltazisq*onetwelfth, &
- deltazisq*4._dp*onethird, &
- -deltazisq*2.5_dp, &
- deltazisq*4._dp*onethird, &
- -deltazisq*onetwelfth /)
-
- f_z(izs) = coef(-2)*f(ize-1) + coef(-1)*f(ize) + coef(0)*f(izs) + coef(1)*f(izs+1) + coef(2)*f(izs+2)
- f_z(izs+1) = coef(-2)*f(ize) + coef(-1)*f(izs) + coef(0)*f(izs+1) + coef(1)*f(izs+2) + coef(2)*f(izs+3)
- do iz=izs+2,ize-2
- f_z(iz) = coef(-2)*f(iz-2) + coef(-1)*f(iz-1) + coef(0)*f(iz) + coef(1)*f(iz+1) + coef(2)*f(iz+2)
- end do
- f_z(ize-1) = coef(-2)*f(ize-3) + coef(-1)*f(ize-2) + coef(0)*f(ize-1) + coef(1)*f(ize) + coef(2)*f(izs)
- f_z(ize) = coef(-2)*f(ize-2) + coef(-1)*f(ize-1) + coef(0)*f(ize) + coef(1)*f(izs) + coef(2)*f(izs+1)
-
-END SUBROUTINE gradzz_fd4
-
-
-
-! Interpolation
-
-SUBROUTINE interp_n2v_fd4( f_n , f_v ) ! implemented for periodic boundary conditions
-
- use prec_const
- use space_grid
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f_n
- real(dp), dimension(izs:ize),intent(out) :: f_v
- real(dp), dimension(-2:1) :: coef ! fd4 stencil
- integer :: iz
-
- coef = (/ -onesixteenth, &
- ninesixteenths, &
- ninesixteenths, &
- -onesixteenth /)
-
- f_v(izs) = coef(-2)*f_n(ize-1) + coef(-1)*f_n(ize) + coef(0)*f_n(izs) + coef(1)*f_n(izs+1)
- f_v(izs+1) = coef(-2)*f_n(ize) + coef(-1)*f_n(izs) + coef(0)*f_n(izs+1) + coef(1)*f_n(izs+2)
- do iz=izs+2,ize-1
- f_v(iz) = coef(-2)*f_n(iz-2) + coef(-1)*f_n(iz-1) + coef(0)*f_n(iz) + coef(1)*f_n(iz+1)
- end do
- f_v(ize) = coef(-2)*f_n(ize-2) + coef(-1)*f_n(ize-1) + coef(0)*f_n(ize) + coef(1)*f_n(izs)
-
-END SUBROUTINE interp_n2v_fd4
-
-
-
-SUBROUTINE interp_v2n_fd4( f_v , f_n ) ! implemented for periodic boundary conditions
-
- use prec_const
- use space_grid
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f_v
- real(dp), dimension(izs:ize),intent(out) :: f_n
- real(dp), dimension(-1:2) :: coef ! fd4 stencil
- integer :: iz
-
- coef = (/ -onesixteenth, &
- ninesixteenths, &
- ninesixteenths, &
- -onesixteenth /)
-
- f_n(izs) = coef(-1)*f_v(ize) + coef(0)*f_v(izs) + coef(1)*f_v(izs+1) + coef(2)*f_v(izs+2)
- do iz=izs+1,ize-2
- f_n(iz) = coef(-1)*f_v(iz-1) + coef(0)*f_v(iz) + coef(1)*f_v(iz+1) + coef(2)*f_v(iz+2)
- end do
- f_n(ize-1) = coef(-1)*f_v(ize-2) + coef(0)*f_v(ize-1) + coef(1)*f_v(ize) + coef(2)*f_v(izs)
- f_n(ize) = coef(-1)*f_v(ize-1) + coef(0)*f_v(ize) + coef(1)*f_v(izs) + coef(2)*f_v(izs+1)
-
-END SUBROUTINE interp_v2n_fd4
-
-
diff --git a/src/grad_up2.h b/src/grad_up2.h
deleted file mode 100644
index 84932368e3b0d8bd2e1d052529926d5f84fb20af..0000000000000000000000000000000000000000
--- a/src/grad_up2.h
+++ /dev/null
@@ -1,33 +0,0 @@
-! For the upwind scheme, the gradients using the two possible sides
-SUBROUTINE gradz_up2_plus( f , f_z )
-
- use prec_const
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f
- real(dp), dimension(izs:ize),intent(out) :: f_z
- integer :: iz
-
- do iz=izs,ize-2
- f_z(iz) = deltazih * ( -3._dp*f(iz) + 4._dp*f(iz+1) -f(iz+2) )
- end do
- f_z(ize-1) = deltazih * ( -3._dp*f(ize-1) + 4._dp*f(ize) -f(izs) )
- f_z(ize) = deltazih * ( -3._dp*f(ize) + 4._dp*f(izs) -f(izs+1) )
-
-END SUBROUTINE gradz_up2_plus
-
-
-SUBROUTINE gradz_up2_minus( f , f_z )
-
- use prec_const
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f
- real(dp), dimension(izs:ize),intent(out) :: f_z
- integer :: iz
-
- f_z(izs) = deltazih * ( f(ize-1) - 4._dp*f(ize) + 3._dp*f(izs) )
- f_z(izs+1) = deltazih * ( f(ize) - 4._dp*f(izs) + 3._dp*f(izs+1) )
- do iz=izs+2,ize
- f_z(iz) = deltazih * ( f(iz-2) - 4._dp*f(iz-1) + 3._dp*f(iz) )
- end do
-
-END SUBROUTINE gradz_up2_minus
diff --git a/src/grad_we4.h b/src/grad_we4.h
deleted file mode 100644
index 97f246cf0ce82f8b40ab91c43612d6540848b504..0000000000000000000000000000000000000000
--- a/src/grad_we4.h
+++ /dev/null
@@ -1,191 +0,0 @@
-! see "Efficient Implementation of Weighted ENO Schemes", GUANG-SHAN JIANG et all, 1996
-
-SUBROUTINE gradz_we4( f , f_z ) ! implemented for periodic boundary conditions
-
- use prec_const
- use space_grid
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f
- real(dp), dimension(izs:ize),intent(out) :: f_z
- real(dp) :: tmp
- integer :: i
-
- call compute_f_hat_plus_one_half_we4( f, f_z )
- ! f_z contains f_hat_plus_one_half
-
- tmp = f_z(ize)
-
- do i=1,ize-izs
- f_z(ize+1-i) = ( f_z(ize+1-i) - f_z(ize-i) )*deltazi
- enddo
- f_z(izs) = ( f_z(izs)- tmp )*deltazi
-
-!! write(*,*) "gradz_we4"
-
-END SUBROUTINE gradz_we4
-
-
-
-SUBROUTINE compute_IS_we4( f, IS )
-
- use prec_const
- use space_grid
- use array
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f
- real(dp), dimension(izs:ize,0:2),intent(out) :: IS
- integer :: iz
-
-
- ! SI_0
- tmp1(izs) = f(ize-1) - 2.0_dp*f(ize) + f(izs)
- tmp1(izs+1) = f(ize) - 2.0_dp*f(izs) + f(izs+1)
- do iz=izs+2,ize
- tmp1(iz) = f(iz-2) - 2.0_dp*f(iz-1) + f(iz)
- enddo
-
- tmp2(izs) = f(ize-1) - 4.0_dp*f(ize) + 3.0_dp*f(izs)
- tmp2(izs+1) = f(ize) - 4.0_dp*f(izs) + 3.0_dp*f(izs+1)
- do iz=izs+2,ize
- tmp2(iz) = f(iz-2) - 4.0_dp*f(iz-1) + 3.0_dp*f(iz)
- enddo
-
- do iz=izs,ize
- IS(iz,0) = thirteentwelfths*tmp1(iz)*tmp1(iz) + 0.25_dp*tmp2(iz)*tmp2(iz)
- enddo
-
-
- ! SI_1
- tmp1(izs) = f(ize) - 2.0_dp*f(izs) + f(izs+1)
- do iz=izs+1,ize-1
- tmp1(iz) = f(iz-1) - 2.0_dp*f(iz) + f(iz+1)
- enddo
- tmp1(ize) = f(ize-1) - 2.0_dp*f(ize) + f(izs)
-
- tmp2(izs) = f(ize) - f(izs+1)
- do iz=izs+1,ize-1
- tmp2(iz) = f(iz-1) - f(iz+1)
- enddo
- tmp2(ize) = f(ize-1)- f(izs)
-
- do iz=izs,ize
- IS(iz,1) = thirteentwelfths*tmp1(iz)*tmp1(iz) + 0.25_dp*tmp2(iz)*tmp2(iz)
- enddo
-
-
- ! SI_2
- do iz=izs,ize-2
- tmp1(iz) = f(iz) - 2.0_dp*f(iz+1) + f(iz+2)
- enddo
- tmp1(ize-1) = f(ize-1) - 2.0_dp*f(ize) + f(izs)
- tmp1(ize) = f(ize) - 2.0_dp*f(izs) + f(izs+1)
-
- do iz=izs,ize-2
- tmp2(iz) = 3.0_dp*f(iz) - 4.0_dp*f(iz+1) + f(iz+2)
- enddo
- tmp2(ize-1) = 3.0_dp*f(ize-1) - 4.0_dp*f(ize) + f(izs)
- tmp2(ize) = 3.0_dp*f(ize) - 4.0_dp*f(izs) + f(izs+1)
-
- do iz=izs,ize
- IS(iz,2) = thirteentwelfths*tmp1(iz)*tmp1(iz) + 0.25_dp*tmp2(iz)*tmp2(iz)
- enddo
-
-END SUBROUTINE compute_IS_we4
-
-
-SUBROUTINE compute_omegas_we4( f, omegas )
-
- use prec_const
- use space_grid
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f
- real(dp), dimension(izs:ize,0:2),intent(out) :: omegas
- real(dp) :: epsi = 0.000001_dp ! 1e-6
- ! real(dp) :: p = 2.0_dp ! this is hardcoded to go faster
- real(dp), dimension(0:2) :: C3k = (/ 0.1, 0.6, 0.3 /)
- integer :: iz
- integer :: k
- real(dp) :: sumalphas ! tmp variable
-
-
- call compute_IS_we4( f, omegas )
- ! omegas now contains IS
-
-
- do k=0,2
- do iz=izs,ize
- omegas(iz,k) = C3k(k) / ( (epsi + omegas(iz,k))*(epsi + omegas(iz,k)) )
-! omegas(iz,k) = C3k(k) / power( epsi + IS, p)
- enddo
- enddo
- ! omegas now contains alphas
-
-
- do iz=izs,ize
- sumalphas = omegas(iz,0)+omegas(iz,1)+omegas(iz,2)
-
- do k=0,2
- omegas(iz,k) = omegas(iz,k) / sumalphas
- enddo
- enddo
- ! omegas now contains omegas
-
-END SUBROUTINE compute_omegas_we4
-
-
-
-
-SUBROUTINE compute_f_hat_plus_one_half_we4( f, f_hat )
-
- use prec_const
- use space_grid
- use array
- implicit none
- real(dp), dimension(izs:ize),intent(in) :: f
- real(dp), dimension(izs:ize),intent(out) :: f_hat
- real(dp) :: epsi = 0.000001_dp ! 1e-6
- ! real(dp) :: p = 2.0_dp ! this is hardcoded to go faster
- real(dp), dimension(0:2,0:2) :: a3kl = reshape((/ onethird, -sevensixths, elevensixths, &
- -onesixth, fivesixths, onethird, &
- onethird, fivesixths, -onesixth /) &
- , shape(a3kl), order=(/2,1/) )
- integer :: iz
-
-
- call compute_omegas_we4( f, omegak_times_q3k )
- ! omegak_times_q3k contains omegas
-
-
- ! k=0
- omegak_times_q3k(izs,0) = omegak_times_q3k(izs,0) *( a3kl(0,0)*f(ize-1)+ a3kl(0,1)*f(ize) + a3kl(0,2)*f(izs) )
- omegak_times_q3k(izs+1,0)= omegak_times_q3k(izs+1,0)*( a3kl(0,0)*f(ize) + a3kl(0,1)*f(izs) + a3kl(0,2)*f(izs+1) )
- do iz=izs+2,ize
- omegak_times_q3k(iz,0) = omegak_times_q3k(iz,0) *( a3kl(0,0)*f(iz-2) + a3kl(0,1)*f(iz-1) + a3kl(0,2)*f(iz) )
- enddo
-
-
- ! k=1
- omegak_times_q3k(izs,1) = omegak_times_q3k(izs,1) *( a3kl(1,0)*f(ize) + a3kl(1,1)*f(izs) + a3kl(1,2)*f(izs+1) )
- do iz=izs+1,ize-1
- omegak_times_q3k(iz,1) = omegak_times_q3k(iz,1) *( a3kl(1,0)*f(iz-1) + a3kl(1,1)*f(iz) + a3kl(1,2)*f(iz+1) )
- enddo
- omegak_times_q3k(ize,1) = omegak_times_q3k(ize,1) *( a3kl(1,0)*f(ize-1)+ a3kl(1,1)*f(ize) + a3kl(1,2)*f(izs) )
-
- ! k=2
- do iz=izs,ize-2
- omegak_times_q3k(iz,2) = omegak_times_q3k(iz,2) *( a3kl(2,0)*f(iz) + a3kl(2,1)*f(iz+1) + a3kl(2,2)*f(iz+2) )
- enddo
- omegak_times_q3k(ize-1,2)= omegak_times_q3k(ize-1,2)*( a3kl(2,0)*f(ize-1)+ a3kl(2,1)*f(ize) + a3kl(2,2)*f(izs) )
- omegak_times_q3k(ize,2) = omegak_times_q3k(ize,2) *( a3kl(2,0)*f(ize) + a3kl(2,1)*f(izs) + a3kl(2,2)*f(izs+1) )
-
-
-
- do iz=izs,ize
- f_hat(iz) = omegak_times_q3k(iz,0) + omegak_times_q3k(iz,1) + omegak_times_q3k(iz,2)
- enddo
-
-END SUBROUTINE compute_f_hat_plus_one_half_we4
-
-
-
-
diff --git a/src/gradients.F90 b/src/gradients.F90
deleted file mode 100644
index a939697eea6c6219a2ca72a08e106bc2d98d92ae..0000000000000000000000000000000000000000
--- a/src/gradients.F90
+++ /dev/null
@@ -1,103 +0,0 @@
-module gradients
- use space_grid
- use prec_const
- implicit none
-
- ! Staggered grids are used, the v-grid of vpar is shifted by deltaz/2 compared to the normal n-grid of theta, temp, moments, phi.
- ! (z_i) on n-grid is (z_i+deltaz/2) of v-grid.
- ! When evaluating a gradient on another grid, we use formulas that interpolate directly.
-
- ! Gradients
- procedure(gradsub), pointer :: gradz => null()
- procedure(gradsub), pointer :: gradz_n2v => null()
- procedure(gradsub), pointer :: gradz_v2n => null()
-
- procedure(gradsub), pointer :: gradzm => null() ! For upwind scheme, other set of spatial derivatives, the "minus side"
-
- procedure(gradsub), pointer :: gradzz => null()
-
-
- ! Grid-grid interpolations
- procedure(gradsub), pointer :: interp_n2v => null()
- procedure(gradsub), pointer :: interp_v2n => null()
-
-
-
- !Interface for gradient subroutines
- abstract interface
- subroutine gradsub( f , fp )
- use space_grid, only: izs,ize
- use prec_const
- implicit none
-
- real(dp), dimension(izs:ize), intent(in) :: f
- real(dp), dimension(izs:ize), intent(out) :: fp
-
- end subroutine gradsub
- end interface
-
-
-contains
- !To promote cleanliness, use .h include files
- !Please put things where they belong if you add new routines
-
- !Inside this file place formulas that are order independent,
- !e.g. sums of gradients and such
-
- !Note that eventually all the functions here should become private
-
-#include "grad_fa2.h"
-#include "grad_up2.h"
-#include "grad_fd4.h"
-#include "grad_we4.h"
-
-
- !In this subroutine, which should be called before any gradient is carried out,
- !we select the numerical scheme using function pointers
- subroutine set_gradient_scheme(gradient_scheme)
-
- use prec_const
- implicit none
- character(len=3),intent(in) :: gradient_scheme
-
- select case ( gradient_scheme )
-! NOT IMPLEMENTED YET, need finite difference order 2 formulas for staggered grids
-! case ( 'fa2' )
-! gradz => gradz_fa2
-! gradz_n2v => gradz_n2v_fa2
-! gradz_v2n => gradz_v2n_fa2
-! gradzz => gradzz_fa2
-! case ( 'up2' )
-! gradz => gradz_up2_plus
-! gradzm => gradz_up2_minus
- case ('fd4')
- ! Gradients
-! gradz => gradz_forfd4
- gradz => gradz_fd4
- gradz_n2v => gradz_n2v_fd4
- gradz_v2n => gradz_v2n_fd4
- gradzz => gradzz_fd4
-
- ! Interpolation
- interp_n2v => interp_n2v_fd4
- interp_v2n => interp_v2n_fd4
-
- case ('we4') ! WENO weighted essentially non-oscillatory order 4, meaning r=3
- ! Gradients
- gradz => gradz_we4
- gradz_n2v => gradz_n2v_fd4
- gradz_v2n => gradz_v2n_fd4
- gradzz => gradzz_fd4
-
- ! Interpolation
- interp_n2v => interp_n2v_fd4
- interp_v2n => interp_v2n_fd4
-
-
- case default
- write (*,*) "Unknown numerical scheme ", gradient_scheme
- end select
-
- end subroutine set_gradient_scheme
-
-end module gradients
diff --git a/src/inital.F90 b/src/inital.F90
index 1c615a8748356c27042013bf9a15940a4d3ccd58..2699fc6c57e7829878c60f05bc48aefabb80315c 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -16,55 +16,34 @@ SUBROUTINE init_profiles
! Set initial conditions
USE basic
- USE space_grid
+ USE fourier_grid
USE fields
use initial_par
USE time_integration
- USE model, ONLY: gradient_scheme
use prec_const
IMPLICIT NONE
- INTEGER :: ip, ij, ikr, ikz
- integer, dimension(12) :: iseedarr
- real(dp) :: noiser, noisez
+ INTEGER :: ipj, ikr, ikz
+ INTEGER, DIMENSION(12) :: iseedarr
+ REAL(dp) :: noise
! Seed random number generator
iseedarr(:)=iseed
CALL RANDOM_SEED(PUT=iseedarr)
-
CALL set_updatetlevel(1)
-
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
-
- CALL RANDOM_NUMBER(noiser)
- CALL RANDOM_NUMBER(noisez)
- theta(ikr,ikz,:) = log( initback_density &
- + init_ampli_density * sin(init_nb_oscil_density*TWOPI/(krmax-krmin)*rarray(ikr)) + initnoise_density*(noiser-0.5_dp) & ! sine profile with noise
- + init_ampli_density * sin(init_nb_oscil_density*TWOPI/(kzmax-kzmin)*zarray(ikz)) + initnoise_density*(noisez-0.5_dp) ) ! sine profile with noise
- CALL RANDOM_NUMBER(noiser)
- CALL RANDOM_NUMBER(noisez)
- temp(ikr,ikz,:) = log( initback_temp + init_ampli_temp * sin(init_nb_oscil_temp*TWOPI/(zmax-zmin)*zarray(iz)) + initnoise_temp*(noise-0.5_dp) ) ! sine profile with noise
-
- CALL RANDOM_NUMBER(noiser)
- CALL RANDOM_NUMBER(noisez)
- vpar(ikr,ikz,:) = initback_vpar + initnoise_vpar*(noisez-0.5_dp)
-
- DO ip=ips,ipe
- DO ij=ijs,ije
- CALL RANDOM_NUMBER(noiser)
- CALL RANDOM_NUMBER(noisez)
- moments(ip,ij,ikr,ikz,:) = initback_moments + initnoise_moments*(noisez-0.5_dp)
- END DO
+ DO ipj=ipjs,ipje
+ CALL RANDOM_NUMBER(noise)
+ moments( ipj, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
END DO
-
END DO
END DO
- call poisson ! To set phi
+ CALL poisson ! To set phi
END SUBROUTINE init_profiles
diff --git a/src/initial_par_mod.F90 b/src/initial_par_mod.F90
index 564d0654aedc2648306b27052533008ad1183d26..64953baa87bfeeb6165976ab07c054e279294740 100644
--- a/src/initial_par_mod.F90
+++ b/src/initial_par_mod.F90
@@ -6,15 +6,9 @@ MODULE initial_par
PRIVATE
! Initial background level
- REAL(dp), PUBLIC, PROTECTED :: initback_density=0._dp
- REAL(dp), PUBLIC, PROTECTED :: initback_temp=0._dp
- REAL(dp), PUBLIC, PROTECTED :: initback_vpar=0._dp
REAL(dp), PUBLIC, PROTECTED :: initback_moments=0._dp
! Initial background noise amplitude
- REAL(dp), PUBLIC, PROTECTED :: initnoise_density=1E-6_dp
- REAL(dp), PUBLIC, PROTECTED :: initnoise_temp=1E-6_dp
- REAL(dp), PUBLIC, PROTECTED :: initnoise_vpar=1E-6_dp
REAL(dp), PUBLIC, PROTECTED :: initnoise_moments=1E-6_dp
! Initialization for random number generator
@@ -38,10 +32,9 @@ CONTAINS
USE prec_const
IMPLICIT NONE
- NAMELIST /INITIAL_CON/ initback_density, initback_temp, initback_vpar, initback_moments
- NAMELIST /INITIAL_CON/ initnoise_density, initnoise_temp, initnoise_vpar, initnoise_moments
+ NAMELIST /INITIAL_CON/ initback_moments
+ NAMELIST /INITIAL_CON/ initnoise_moments
NAMELIST /INITIAL_CON/ iseed
- NAMELIST /INITIAL_CON/ init_nb_oscil_density, init_nb_oscil_temp, init_ampli_density, init_ampli_temp
READ(lu_in,initial_con)
WRITE(*,initial_con)
@@ -58,23 +51,12 @@ CONTAINS
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
- CALL attach(fidres, TRIM(str), "initback_theta", initback_density)
- CALL attach(fidres, TRIM(str), "initback_temp", initback_temp)
- CALL attach(fidres, TRIM(str), "initback_vpar", initback_vpar)
CALL attach(fidres, TRIM(str), "initback_moments", initback_moments)
- CALL attach(fidres, TRIM(str), "initnoise_theta", initnoise_density)
- CALL attach(fidres, TRIM(str), "initnoise_temp", initnoise_temp)
- CALL attach(fidres, TRIM(str), "initnoise_vpar", initnoise_vpar)
CALL attach(fidres, TRIM(str), "initnoise_moments", initnoise_moments)
CALL attach(fidres, TRIM(str), "iseed", iseed)
- CALL attach(fidres, TRIM(str), "init_nb_oscil_theta", init_nb_oscil_density)
- CALL attach(fidres, TRIM(str), "init_nb_oscil_temp", init_nb_oscil_temp)
- CALL attach(fidres, TRIM(str), "init_ampli_theta", init_ampli_density)
- CALL attach(fidres, TRIM(str), "init_ampli_temp", init_ampli_temp)
-
END SUBROUTINE initial_outputinputs
END MODULE initial_par
diff --git a/src/memory.F90 b/src/memory.F90
index 94ada84a01f92acb2490597adb8d715970a919b8..035697f0528d56ddfc4bb2a0e947b4a31ec5ac67 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -4,8 +4,7 @@ SUBROUTINE memory
USE array
USE basic
USE fields
- USE model, ONLY: gradient_scheme
- USE space_grid
+ USE fourier_grid
USE time_integration
use prec_const
@@ -13,109 +12,20 @@ SUBROUTINE memory
- ! Principal fields
- CALL allocate_array(theta,izs,ize,1,ntimelevel)
- CALL allocate_array(temp,izs,ize,1,ntimelevel)
- CALL allocate_array(vpar,izs,ize,1,ntimelevel)
- CALL allocate_array(moments,ips,ipe,izs,ize,1,ntimelevel)
- CALL allocate_array(phi,izs,ize)
- ! On opposite grid
- CALL allocate_array(vpar_n,izs,ize)
-
- ! Right hand sides of the time evolution equations
- CALL allocate_array(theta_rhs,izs,ize,1,ntimelevel)
- CALL allocate_array(temp_rhs,izs,ize,1,ntimelevel)
- CALL allocate_array(vpar_rhs,izs,ize,1,ntimelevel)
- CALL allocate_array(moments_rhs,ips,ipe,izs,ize,1,ntimelevel)
- ! On opposite grid
- CALL allocate_array(temp_rhs_v,izs,ize)
- CALL allocate_array(vpar_rhs_n,izs,ize)
-
- ! Auxiliary quantities
- CALL allocate_array(sqrt_exp_temp,izs,ize)
- CALL allocate_array(sqrt_exp_temp_v,izs,ize)
-
- ! Spatial 1rst derivatives on respective grids
- CALL allocate_array(thetaz,izs,ize)
- CALL allocate_array(sqrt_exp_tempz,izs,ize)
- CALL allocate_array(vparz,izs,ize)
- CALL allocate_array(phiz,izs,ize)
- CALL allocate_array(momentsz,ips,ipe,izs,ize)
- ! On opposite grid
- CALL allocate_array(thetaz_v,izs,ize)
- CALL allocate_array(sqrt_exp_tempz_v,izs,ize)
- CALL allocate_array(vparz_n,izs,ize)
- CALL allocate_array(phiz_v,izs,ize)
-
-
- ! select case (gradient_scheme)
- ! case('up2') ! If upwind scheme : need spatial derivaties for the "minus side"
- ! CALL allocate_array(thetazm,izs,ize)
- ! CALL allocate_array(tempzm,izs,ize)
- ! CALL allocate_array(vparzm,izs,ize)
- ! CALL allocate_array(phizm,izs,ize)
- ! CALL allocate_array(momentszm,ips,ipe,izs,ize)
-
- ! CALL allocate_array(thetazz,0,0)
- ! CALL allocate_array(tempzz,0,0)
- ! CALL allocate_array(vparzz,0,0)
- ! CALL allocate_array(momentszz,0,0,0,0)
-
- ! case('fa2','fd4') ! If no upwind scheme : numerical diffusion needs 2nd order spatial derivatives
- ! CALL allocate_array(thetazm,0,0)
- ! CALL allocate_array(tempzm,0,0)
- ! CALL allocate_array(vparzm,0,0)
- ! CALL allocate_array(phizm,0,0)
- ! CALL allocate_array(momentszm,0,0,0,0)
-
- CALL allocate_array(thetazz,izs,ize)
- CALL allocate_array(tempzz,izs,ize)
- CALL allocate_array(vparzz,izs,ize)
- CALL allocate_array(momentszz,ips,ipe,izs,ize)
- ! end select
+ ! Moments and moments rhs
+ CALL allocate_array( moments, ipjs,ipje, ikrs,ikre, ikzs,ikze, 1,ntimelevel)
+ CALL allocate_array( moments_rhs, ipjs,ipje, ikrs,ikre, ikzs,ikze, 1,ntimelevel)
+ ! Electrostatic potential
+ CALL allocate_array( phi, ikrs,ikre, ikzs,ikze)
! Intermediate steps in rhs of equations
- ! theta
- CALL allocate_array(theta_Cpl,izs,ize)
- CALL allocate_array(theta_Hpl,izs,ize)
- ! temp
- CALL allocate_array(temp_Cpl,izs,ize)
- CALL allocate_array(temp_Fpl,izs,ize)
- CALL allocate_array(temp_Hpl,izs,ize)
- CALL allocate_array(temp_Ipl,izs,ize)
- ! vpar
- CALL allocate_array(vpar_Cpl,izs,ize)
- CALL allocate_array(vpar_Dpl,izs,ize)
- CALL allocate_array(vpar_Epl,izs,ize)
- CALL allocate_array(vpar_Fpl,izs,ize)
- CALL allocate_array(vpar_Hpl,izs,ize)
- ! moments
- CALL allocate_array(moments_Apl,izs,ize)
- CALL allocate_array(moments_Bpl,izs,ize)
- CALL allocate_array(moments_Cpl,izs,ize)
- CALL allocate_array(moments_Dpl,izs,ize)
- CALL allocate_array(moments_Epl,izs,ize)
- CALL allocate_array(moments_Fpl,izs,ize)
- CALL allocate_array(moments_Gpl,izs,ize)
- CALL allocate_array(moments_Hpl,izs,ize)
-
-
- ! Poisson solver
- ! CALL init(neq, 1, laplace_smumps)
- CALL allocate_array(laplace_rhs,1,neq)
- CALL allocate_array(laplace_sol,1,neq)
-
-
- ! FFT
- CALL allocate_array(fft_filter,izs,ize)
- CALL allocate_array(fft_input,izs,ize)
-
-
- ! WENO gradients
- CALL allocate_array(tmp1,izs,ize)
- CALL allocate_array(tmp2,izs,ize)
- CALL allocate_array(omegak_times_q3k,izs,ize,0,2)
-
-
+ !CALL allocate_array( moments_Apl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
+ !CALL allocate_array( moments_Bpl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
+ !CALL allocate_array( moments_Cpl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
+ !CALL allocate_array( moments_Dpl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
+ !CALL allocate_array( moments_Epl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
+ !CALL allocate_array( moments_Fpl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
+ !CALL allocate_array( moments_Gpl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
+ !CALL allocate_array( moments_Hpl, ipjs,ipje, ikrs,ikre, ikzs,ikze,)
END SUBROUTINE memory
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index e24dbfe50861a702ff2473b55dfefcbd1f476c3a..fe657e64cd89065f8819fb074be4feb42369dbd0 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -5,28 +5,19 @@ MODULE model
IMPLICIT NONE
PRIVATE
- REAL(dp), PUBLIC, PROTECTED :: nu = 1._dp ! Collision frequency
-
- ! (Numerical) diffusion coefficient
- REAL(dp), PUBLIC, PROTECTED :: diff_theta = 1._dp
- REAL(dp), PUBLIC, PROTECTED :: diff_temp = 1._dp
- REAL(dp), PUBLIC, PROTECTED :: diff_vpar = 1._dp
- REAL(dp), PUBLIC, PROTECTED :: diff_moments = 1._dp
-
- ! Fast Fourier Transform to filter out high frequence
- LOGICAL, PUBLIC, PROTECTED :: fft_suppress_high_freq = .false.
- REAL(dp), PUBLIC, PROTECTED :: fft_sigma = 2._dp
-
-
- CHARACTER(len=3), PUBLIC, PROTECTED :: gradient_scheme='fd4'
-
- LOGICAL, PUBLIC, PROTECTED :: freeze_theta = .false.
- LOGICAL, PUBLIC, PROTECTED :: freeze_temp = .false.
- LOGICAL, PUBLIC, PROTECTED :: freeze_vpar = .false.
- LOGICAL, PUBLIC, PROTECTED :: freeze_moments = .false.
- LOGICAL, PUBLIC, PROTECTED :: freeze_phi = .false.
-
- PUBLIC :: model_readinputs, model_outputinputs
+ REAL(dp), PUBLIC, PROTECTED :: nu = 1._dp ! Collision frequency
+ REAL(dp), PUBLIC, PROTECTED :: tau_e = 1._dp ! Collision frequency
+ REAL(dp), PUBLIC, PROTECTED :: tau_i = 1._dp ! Collision frequency
+ REAL(dp), PUBLIC, PROTECTED :: sigma_e = 1._dp ! Collision frequency
+ REAL(dp), PUBLIC, PROTECTED :: sigma_i = 1._dp ! Collision frequency
+ REAL(dp), PUBLIC, PROTECTED :: q_e = 1._dp ! Collision frequency
+ REAL(dp), PUBLIC, PROTECTED :: q_i = 1._dp ! Collision frequency
+ REAL(dp), PUBLIC, PROTECTED :: eta_n = 1._dp ! Collision frequency
+ REAL(dp), PUBLIC, PROTECTED :: eta_T = 1._dp ! Collision frequency
+ REAL(dp), PUBLIC, PROTECTED :: eta_B = 1._dp ! Collision frequency
+ REAL(dp), PUBLIC, PROTECTED :: lambdaD = 1._dp ! Collision frequency
+
+ PUBLIC :: model_readinputs, model_outputinputs
CONTAINS
@@ -37,11 +28,8 @@ CONTAINS
USE prec_const
IMPLICIT NONE
- NAMELIST /MODEL_PAR/ nu
- NAMELIST /MODEL_PAR/ diff_theta, diff_temp, diff_vpar, diff_moments
- NAMELIST /MODEL_PAR/ gradient_scheme
- NAMELIST /MODEL_PAR/ freeze_theta, freeze_temp, freeze_vpar, freeze_moments, freeze_phi
- NAMELIST /MODEL_PAR/ fft_suppress_high_freq, fft_sigma
+ NAMELIST /MODEL_PAR/ nu, tau_e, tau_i, sigma_e, sigma_i, &
+ q_e, q_i, eta_n, eta_T, eta_B, lambdaD
READ(lu_in,model_par)
WRITE(*,model_par)
@@ -58,23 +46,20 @@ CONTAINS
USE futils, ONLY: attach
USE prec_const
IMPLICIT NONE
+
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
-
- CALL attach(fidres, TRIM(str), "nu", nu)
- CALL attach(fidres, TRIM(str), "diff_theta", diff_theta)
- CALL attach(fidres, TRIM(str), "diff_temp", diff_temp)
- CALL attach(fidres, TRIM(str), "diff_vpar", diff_vpar)
- CALL attach(fidres, TRIM(str), "diff_moments", diff_moments)
- CALL attach(fidres, TRIM(str), "gradient_scheme", gradient_scheme)
- CALL attach(fidres, TRIM(str), "freeze_theta", freeze_theta)
- CALL attach(fidres, TRIM(str), "freeze_temp", freeze_temp)
- CALL attach(fidres, TRIM(str), "freeze_vpar", freeze_vpar)
- CALL attach(fidres, TRIM(str), "freeze_moments", freeze_moments)
- CALL attach(fidres, TRIM(str), "freeze_phi", freeze_phi)
- CALL attach(fidres, TRIM(str), "fft_suppress_high_freq", fft_suppress_high_freq)
- CALL attach(fidres, TRIM(str), "fft_sigma", fft_sigma)
-
+ CALL attach(fidres, TRIM(str), "nu", nu)
+ CALL attach(fidres, TRIM(str), "tau_e", tau_e)
+ CALL attach(fidres, TRIM(str), "tau_i", tau_i)
+ CALL attach(fidres, TRIM(str), "sigma_e", sigma_e)
+ CALL attach(fidres, TRIM(str), "sigma_i", sigma_i)
+ CALL attach(fidres, TRIM(str), "q_e", q_e)
+ CALL attach(fidres, TRIM(str), "q_i", q_i)
+ CALL attach(fidres, TRIM(str), "eta_n", eta_n)
+ CALL attach(fidres, TRIM(str), "eta_T", eta_T)
+ CALL attach(fidres, TRIM(str), "eta_B", eta_B)
+ CALL attach(fidres, TRIM(str), "lambdaD", lambdaD)
END SUBROUTINE model_outputinputs
END MODULE model
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 83174fbb0b61fa2161ba92d9eca8b7b44701a06c..86172659b86d597637d074394668e0b5cd6b81ab 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -4,147 +4,48 @@ SUBROUTINE moments_eq_rhs
USE time_integration
USE array
USE fields
- USE space_grid
+ USE fourier_grid
USE model
use prec_const
IMPLICIT NONE
- INTEGER:: ip, iz
- REAL(dp) :: ip_dp
+ INTEGER:: ip,ij, ipj, ikr,ikz
+ REAL(dp) :: ip_dp, ij_dp
REAL(dp) :: tmp, moml, sqrtT
- do ip = ips, ipe
- ip_dp = real(ip,dp) ! From int to double (compute once)
+ do ipj = ipjs, ipje
+ if (ipj .le. Nmome + 1) then ! electrons moments
+ CALL rabe(ipj, ip, ij) ! Compute p,j electrons moments degree
+ else ! ions moments
+ CALL rabi(ipj, ip, ij) ! Compute p,j ions moments degree
+ endif
- do iz = izs,ize ! Compute coefficients
+ ip_dp = real(ip,dp)
+ ij_dp = real(ij,dp)
- sqrtT = sqrt_exp_temp(iz)
+ do ikr = ikrs,ikre
+ do ikz = ikzs,ikze
! N_e^{p} term
- moml = moments(ip,iz,updatetlevel)
- moments_Apl(iz) = -nu*ip_dp*moml ! Collisional damping from Lenard-Bernstein Operator
- moments_Bpl(iz) = -moml
- moments_Cpl(iz) = -sqrtT*vpar_n(iz)*INVSQRT2*moml
- moments_Dpl(iz) = 0._dp
- moments_Epl(iz) = -ip_dp*0.5_dp*moml
- moments_Fpl(iz) = -vpar_n(iz)*2._dp*SQRT2*ip_dp*moml
-! moments_Fpl(iz) = -sqrtT*vpar_n(iz)*SQRT2*ip_dp*moml
- moments_Gpl(iz) = 0._dp
- moments_Hpl(iz) = -sqrtT*SQRT2*(ip_dp+1._dp)*moml
! N_e^{p+1} term
- if (ip+1 .le. ipe) then
- moml = moments(ip+1,iz,updatetlevel)
- moments_Cpl(iz) = moments_Cpl(iz) -sqrtT*sqrt(ip_dp+1._dp)*0.5_dp*moml
- moments_Fpl(iz) = moments_Fpl(iz) -ip_dp*sqrt(ip_dp+1._dp)*moml
-! moments_Fpl(iz) = moments_Fpl(iz) -sqrtT*ip_dp*sqrt(ip_dp+1._dp)*0.5_dp*moml
- endif
! N_e^{p-1} term
- if (ip-1 .ge. ips) then
- moml = moments(ip-1,iz,updatetlevel)
- if (ip .ne. 3) then
- moments_Apl(iz) = moments_Apl(iz) -nu*SQRT2*sqrt(ip_dp)*moml*vpar_n(iz) ! Collisions
- endif
- moments_Cpl(iz) = moments_Cpl(iz) -sqrtT*sqrt(ip_dp)*0.5_dp*moml
- moments_Dpl(iz) = moments_Dpl(iz) +sqrt(ip_dp)/sqrtT*moml
- moments_Epl(iz) = moments_Epl(iz) -vpar_n(iz)*sqrt(ip_dp*0.5_dp)*moml
- moments_Fpl(iz) = moments_Fpl(iz) -sqrt(ip_dp)*(2._dp*ip_dp-1._dp+2._dp*vpar_n(iz)*vpar_n(iz))*moml
-! moments_Fpl(iz) = moments_Fpl(iz) -sqrtT*sqrt(ip_dp)*0.5_dp*(2._dp*ip_dp-1._dp+2._dp*vpar_n(iz)*vpar_n(iz))*moml
- moments_Gpl(iz) = moments_Gpl(iz) -SQRT2*sqrt(ip_dp)*moml
- moments_Hpl(iz) = moments_Hpl(iz) -sqrtT*vpar_n(iz)*2._dp*sqrt(ip_dp)*moml
- endif
! N_e^{p-2} term
- if (ip-2 .ge. ips) then
- moml = moments(ip-2,iz,updatetlevel)
- moments_Epl(iz) = moments_Epl(iz) -sqrt(ip_dp*(ip_dp-1._dp))*0.5_dp*moml
- moments_Fpl(iz) = moments_Fpl(iz) -2._dp*sqrt(2._dp*ip_dp*(ip_dp-1._dp))*moml
-! moments_Fpl(iz) = moments_Fpl(iz) -sqrtT*sqrt(2._dp*ip_dp*(ip_dp-1._dp))*moml
- moments_Hpl(iz) = moments_Hpl(iz) -sqrtT*sqrt(2._dp*ip_dp*(ip_dp-1._dp))*moml
- endif
-
-
- ! N_e^{p-3} term
- if (ip-3 .ge. ips) then
- moments_Fpl(iz) = moments_Fpl(iz) -sqrt(ip_dp*(ip_dp-1._dp)*(ip_dp-2._dp))*moments(ip-3,iz,updatetlevel)
-! moments_Fpl(iz) = moments_Fpl(iz) -sqrtT*sqrt(ip_dp*(ip_dp-1._dp)*(ip_dp-2._dp))*0.5_dp*moments(ip-3,iz,updatetlevel)
- else if (ip .eq. ips) then ! ip-3=0 so N_e^0 is equal to 1
- moments_Fpl(iz) = moments_Fpl(iz) -SQRT3*SQRT2
-! moments_Fpl(iz) = moments_Fpl(iz) -sqrtT*SQRT3*SQRT2*0.5_dp
- endif
end do
-
- select case (gradient_scheme)
- ! case('up2') ! Upwind scheme, check the sign of coefficients
- ! do iz = izs,ize
- ! if (moments_Cpl(iz)>0) then
- ! tmp = moments_Cpl(iz)*thetaz(iz)
- ! else
- ! tmp = moments_Cpl(iz)*thetazm(iz)
- ! endif
-
- ! if (moments_Dpl(iz)>0) then
- ! tmp = tmp + moments_Dpl(iz)*phiz(iz)
- ! else
- ! tmp = tmp + moments_Dpl(iz)*phizm(iz)
- ! endif
-
- ! if (moments_Fpl(iz)>0) then
- ! tmp = tmp + moments_Fpl(iz)*tempz(iz)
- ! else
- ! tmp = tmp + moments_Fpl(iz)*tempzm(iz)
- ! endif
-
- ! if (moments_Hpl(iz)>0) then
- ! tmp = tmp + moments_Hpl(iz)*vparz(iz)
- ! else
- ! tmp = tmp + moments_Hpl(iz)*vparzm(iz)
- ! endif
-
- ! if (vpar(iz,updatetlevel)<0) then ! I_p^l
- ! moments_Ipl = -sqrtT*SQRT2*vpar(iz,updatetlevel)*momentsz(ip,iz)
- ! else
- ! moments_Ipl = -sqrtT*SQRT2*vpar(iz,updatetlevel)*momentszm(ip,iz)
- ! endif
- ! if (ip+1 .le. ipe) moments_Ipl = moments_Ipl -sqrtT*sqrt(ip_dp+1._dp)*momentszm(ip+1,iz) ! coefficient always negative
- ! if (ip-1 .ge. ips) moments_Ipl = moments_Ipl -sqrtT*sqrt(ip_dp)*momentszm(ip-1,iz) ! coefficient always negative
-
-
- ! moments_rhs(ip,iz,updatetlevel) = tmp &
- ! +moments_Apl(iz) &
- ! +moments_Bpl(iz)*theta_rhs(iz,updatetlevel) &
- ! +moments_Epl(iz)*temp_rhs(iz,updatetlevel) &
- ! +moments_Gpl(iz)*vpar_rhs(iz,updatetlevel) &
- ! +moments_Ipl
-
- ! end do
- case('fa2','fd4','we4')
- do iz = izs,ize
-
- moments_Ipl = -sqrtT*SQRT2*vpar_n(iz)*momentsz(ip,iz)
- if (ip+1 .le. ipe) moments_Ipl = moments_Ipl -sqrtT*sqrt(ip_dp+1._dp)*momentsz(ip+1,iz) ! coefficient always negative
- if (ip-1 .ge. ips) moments_Ipl = moments_Ipl -sqrtT*sqrt(ip_dp)*momentsz(ip-1,iz) ! coefficient always negative
-
- moments_rhs(ip,iz,updatetlevel) = diff_moments*momentszz(ip,iz) & ! numerical diffusion added for stability
- +moments_Apl(iz) &
- +moments_Bpl(iz)*theta_rhs(iz,updatetlevel) &
- +moments_Cpl(iz)*thetaz(iz) &
- +moments_Dpl(iz)*phiz(iz) &
- +moments_Epl(iz)*temp_rhs(iz,updatetlevel) &
- +moments_Fpl(iz)*sqrt_exp_tempz(iz) &
- +moments_Gpl(iz)*vpar_rhs_n(iz) &
- +moments_Hpl(iz)*vparz_n(iz) &
- +moments_Ipl
-! moments_rhs(ip,iz,updatetlevel) = 0._dp ! debug
+ end do
+
+ do ikr = ikrs,ikre
+ do ikz = ikzs,ikze
+ moments_rhs(ipj,ikr,ikz,updatetlevel) = 0._dp ! debug
end do
- end select
-
+ end do
end do
END SUBROUTINE moments_eq_rhs
diff --git a/src/numerics.F90 b/src/numerics.F90
deleted file mode 100644
index e384cde7c02c2475a7e84a92110b5b2478fe09d3..0000000000000000000000000000000000000000
--- a/src/numerics.F90
+++ /dev/null
@@ -1,96 +0,0 @@
-SUBROUTINE evaluation_auxfield_total
-
- USE space_grid
- USE fields
- USE array
- USE time_integration
- USE gradients
- USE prec_const
-
- IMPLICIT NONE
- integer :: iz
- integer :: ip
-
- DO iz = izs, ize
- ! Actual temperature = exp ( ln temp )
- ! Actual density = exp ( ln theta )
- sqrt_exp_temp(iz) = exp(temp(iz,updatetlevel)*0.5_dp)
- END DO
-
- CALL parallel_gradients
-
- CALL interpolation
-
-END SUBROUTINE evaluation_auxfield_total
-
-
-!Compute parallel gradients for RHS of equations
-SUBROUTINE parallel_gradients
- USE array
- USE fields
- USE gradients
- USE model, only: gradient_scheme
- USE space_grid
- USE time_integration, only: updatetlevel
-
- use prec_const
- IMPLICIT NONE
-
- integer :: ip
-
- CALL gradz(theta(:,updatetlevel), thetaz)
- CALL gradz(sqrt_exp_temp, sqrt_exp_tempz)
- CALL gradz(vpar(:,updatetlevel), vparz)
- CALL gradz(phi(:), phiz)
- do ip=ips,ipe
- CALL gradz(moments(ip,:,updatetlevel), momentsz(ip,:))
- enddo
-
- ! Derivatives that are on the opposite grid
- CALL gradz_n2v(theta(:,updatetlevel), thetaz_v)
- CALL gradz_n2v(sqrt_exp_temp, sqrt_exp_tempz_v)
- CALL gradz_v2n(vpar(:,updatetlevel), vparz_n)
- CALL gradz_n2v(phi, phiz_v)
-
- select case (gradient_scheme)
-! case('up2') ! If upwind scheme : need spatial derivaties for the "minus side"
-! CALL gradzm(theta(:,updatetlevel), thetazm)
-! CALL gradzm(temp(:,updatetlevel), tempzm)
-! CALL gradzm(vpar(:,updatetlevel), vparzm)
-! CALL gradzm(phi(:), phizm)
-! do ip=ips,ipe
-! CALL gradzm(moments(ip,:,updatetlevel), momentszm(ip,:))
-! enddo
-
- case('fa2','fd4','we4') ! If no upwind scheme : numerical diffusion needs 2nd order spatial derivatives
- CALL gradzz(theta(:,updatetlevel), thetazz)
- CALL gradzz(temp(:,updatetlevel), tempzz)
- CALL gradzz(vpar(:,updatetlevel), vparzz)
- do ip=ips,ipe
- CALL gradzz(moments(ip,:,updatetlevel), momentszz(ip,:))
- enddo
- end select
-
-
-END SUBROUTINE parallel_gradients
-
-
-
-SUBROUTINE interpolation
- USE array
- USE fields
- USE gradients
- USE time_integration, only: updatetlevel
-
- use prec_const
- IMPLICIT NONE
-
-
- CALL interp_n2v(sqrt_exp_temp, sqrt_exp_temp_v)
-
- CALL interp_v2n(vpar(:,updatetlevel), vpar_n)
-
-
-
-END SUBROUTINE interpolation
-
diff --git a/src/readinputs.F90 b/src/readinputs.F90
index fecaeb2b5ac9a7d2e7752e5b93cfe3576f7666b4..cdcb0cc472a3fe359e28ce8fc45dca032d8118d5 100644
--- a/src/readinputs.F90
+++ b/src/readinputs.F90
@@ -1,7 +1,7 @@
SUBROUTINE readinputs
! Additional data specific for a new run
- USE space_grid, ONLY: space_grid_readinputs
+ USE fourier_grid, ONLY: fourier_grid_readinputs
USE diagnostics_par, ONLY: output_par_readinputs
USE model, ONLY: model_readinputs
USE initial_par, ONLY: initial_readinputs
@@ -15,7 +15,7 @@ SUBROUTINE readinputs
! The input file must be in the same order as the following read routines commands (eg spacegrid then output then model)
! Load grid data from input file
- CALL space_grid_readinputs
+ CALL fourier_grid_readinputs
! Load diagnostic options from input file
CALL output_par_readinputs
diff --git a/src/space_grid_mod.F90 b/src/space_grid_mod.F90
deleted file mode 100644
index 7286c92eff75237bd16a93c3ef74963496103ddb..0000000000000000000000000000000000000000
--- a/src/space_grid_mod.F90
+++ /dev/null
@@ -1,117 +0,0 @@
-MODULE space_grid
- ! Grid module for spatial discretization
-
- USE prec_const
- IMPLICIT NONE
- PRIVATE
-
- ! GRID Namelist
- INTEGER, PUBLIC, PROTECTED :: pmax=3 ! The maximal Hermite-moment computed
- INTEGER, PUBLIC, PROTECTED :: nz=1 ! Number of total internal grid points in z
- REAL(dp), PUBLIC, PROTECTED :: zmin=0._dp ! z coordinate for left boundary
- REAL(dp), PUBLIC, PROTECTED :: zmax=1._dp ! z coordinate for right boundary
-
-
- ! Indices of s -> start , e-> end
- INTEGER, PUBLIC, PROTECTED :: ips, ipe
- INTEGER, PUBLIC, PROTECTED :: izs, ize
-
- ! Toroidal direction
- REAL(dp), PUBLIC, PROTECTED :: deltaz
- real(dp), PUBLIC, PROTECTED :: deltazi
- real(dp), PUBLIC, PROTECTED :: deltazih
- real(dp), PUBLIC, PROTECTED :: deltazisq
-
- ! Grids containing position
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: parray
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: zarray
-
- INTEGER, PUBLIC :: neq ! Number of equations in poisson solver (number of rows/cols of laplace matrix)
-
- ! Public Functions
- PUBLIC :: set_pgrid, set_zgrid
- PUBLIC :: space_grid_readinputs, space_grid_outputinputs
-
-contains
-
- subroutine set_pgrid
- ! Initialize p grid (Hermite moment) parameters
-
- use prec_const
- implicit none
- integer :: ip
-
- ips=3 ! Only start at 3 since N_e^0=1 and N_e^1 = N_e^2 = 0
- ipe=pmax
-
- allocate(parray(ips:ipe))
- DO ip = ips,ipe
- parray(ip) = ip ! simply contains the moment number
- END DO
-
- end subroutine set_pgrid
-
-
- subroutine set_zgrid
- ! Initialize z grid
-
- use prec_const
- implicit none
-
- integer :: iz
-
- ! Start and end indices of grid
- izs=1
- ize=nz
-
- ! Grid spacings, precompute some inverses
- deltaz = (zmax-zmin)/real(nz,dp) !zmax*2._dp*pi/real(nz,dp)
- deltazi = 1._dp/deltaz ! inverse
- deltazih = 0.5_dp/deltaz ! half of inverse
- deltazisq = 1._dp/deltaz/deltaz ! inverse squared
-
- ! Discretized z positions
- allocate(zarray(izs:ize))
- DO iz = izs,ize
- zarray(iz) = zmin + real(iz-1,dp)*deltaz
- END DO
-
- end subroutine set_zgrid
-
-
- SUBROUTINE space_grid_readinputs
- ! Read the input parameters
-
- USE basic, ONLY : lu_in
-
- USE prec_const
- IMPLICIT NONE
-
- NAMELIST /GRID/ pmax, nz, zmin, zmax
-
- READ(lu_in,grid)
- WRITE(*,grid)
-
- END SUBROUTINE space_grid_readinputs
-
-
- SUBROUTINE space_grid_outputinputs(fidres, str)
- ! Write the input parameters to the results_xx.h5 file
-
- USE futils, ONLY: attach
-
- USE prec_const
- IMPLICIT NONE
-
- INTEGER, INTENT(in) :: fidres
- CHARACTER(len=256), INTENT(in) :: str
-
- CALL attach(fidres, TRIM(str), "pmax", pmax)
- CALL attach(fidres, TRIM(str), "nz", nz)
- CALL attach(fidres, TRIM(str), "zmin", zmin)
- CALL attach(fidres, TRIM(str), "zmax", zmax)
-
- END SUBROUTINE space_grid_outputinputs
-
-
-END MODULE space_grid
diff --git a/src/stepon.F90 b/src/stepon.F90
index 9be9f5c58017c88be9412f7a3c4a416d019f692d..db13a317a8cc16167e0252d7c2dbe7c7a9149f24 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -3,65 +3,43 @@ SUBROUTINE stepon
USE basic
USE time_integration
- USE fields, ONLY: theta, temp, vpar, moments, phi
- USE array , ONLY: theta_rhs, temp_rhs, vpar_rhs, moments_rhs
- USE space_grid
+ USE fields, ONLY: moments, phi
+ USE array , ONLY: moments_rhs
+ USE fourier_grid
USE advance_field_routine, ONLY: advance_time_level, advance_field
- USE fft
USE model
USE utility, ONLY: checkfield
use prec_const
IMPLICIT NONE
- INTEGER :: num_step, ip
+ INTEGER :: num_step, ipj
DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
- CALL evaluation_auxfield_total ! Compute gradients, etc.
-
- ! Compute right hand side of derivative the fields
- if (.not. freeze_theta) CALL theta_eq_rhs
- if (.not. freeze_temp) CALL temp_eq_rhs
- if (.not. freeze_vpar) CALL vpar_eq_rhs
- if (.not. freeze_moments) CALL moments_eq_rhs
+ ! Compute right hand side
+ CALL moments_eq_rhs
CALL advance_time_level ! Advance from updatetlevel to updatetlevel+1
- if (.not. freeze_theta) CALL advance_field(theta,theta_rhs)
- if (.not. freeze_temp) CALL advance_field(temp,temp_rhs)
- if (.not. freeze_vpar) CALL advance_field(vpar,vpar_rhs)
- do ip=ips,ipe
- do ij=ijs,ije
- if (.not. freeze_moments) CALL advance_field(moments(ip,ij,:,:,:),moments_rhs(ip,ij,:,:,:))
- enddo
+ do ipj=ipjs,ipje
+ CALL advance_field(moments(ipj,:,:,:),moments_rhs(ipj,:,:,:))
enddo
- if (.not. freeze_phi) CALL poisson
+ ! Solving Poisson equation
+ CALL poisson
CALL checkfield_all()
- END DO
-
- if ( fft_suppress_high_freq ) then ! should have updatetlevel = 1, only suppress high frequences after the runge kutta step
- call suppress_high_freq(theta(:,updatetlevel))
- call suppress_high_freq(temp(:,updatetlevel))
- call suppress_high_freq(vpar(:,updatetlevel))
- do ip=ips,ipe
- call suppress_high_freq(moments(ip,:,updatetlevel))
- enddo
- endif
+ END DO
CONTAINS
SUBROUTINE checkfield_all ! Check all the fields for inf or nan
IF(.NOT.nlend) THEN
- nlend=nlend .or.checkfield(theta(:,updatetlevel),' theta')
- nlend=nlend .or. checkfield(temp(:,updatetlevel),' temp')
- nlend=nlend .or. checkfield(vpar(:,updatetlevel),' vpar')
nlend=nlend .or. checkfield(phi,' phi')
- do ip=ips,ipe
- nlend=nlend .or. checkfield(moments(ip,:,updatetlevel),' moments')
+ do ipj=ipjs,ipje
+ nlend=nlend .or. checkfield(moments(ipj,:,:,updatetlevel),' moments')
enddo
ENDIF
END SUBROUTINE checkfield_all
diff --git a/src/temp_eq_rhs.F90 b/src/temp_eq_rhs.F90
deleted file mode 100644
index 17b21d569689edfa8c03674a2ab3f4c15668c044..0000000000000000000000000000000000000000
--- a/src/temp_eq_rhs.F90
+++ /dev/null
@@ -1,64 +0,0 @@
-SUBROUTINE temp_eq_rhs
-
- USE basic
- USE time_integration
- USE array
- USE fields
- USE space_grid
- USE model
- USE gradients
-
- use prec_const
- IMPLICIT NONE
-
- INTEGER :: iz
- real(dp) :: tmp
-
-
- do iz = izs,ize ! Compute coefficients
- temp_Cpl(iz) = -sqrt_exp_temp(iz)*SQRT3*INVSQRT2*moments(3,iz,updatetlevel)
- temp_Fpl(iz) = -2._dp*SQRT2*(SQRT3*moments(3,iz,updatetlevel)+2._dp)
-! temp_Fpl(iz) = -sqrt_exp_temp(iz)*SQRT2*(SQRT3*moments(3,iz,updatetlevel)+2._dp)
-! temp_Fpl(iz) = -SQRT2*(SQRT3*moments(3,iz,updatetlevel)+2._dp)
- temp_Hpl(iz) = -sqrt_exp_temp(iz)*2._dp*SQRT2
- temp_Ipl(iz) = -sqrt_exp_temp(iz)*SQRT2*SQRT3
- end do
-
-
-
- select case (gradient_scheme)
- ! case('up2') ! Upwind scheme, check the sign of coefficients
- ! do iz = izs,ize
- ! if (temp_Cpl(iz)>0) then
- ! tmp = temp_Cpl(iz)*thetaz(iz)
- ! else
- ! tmp = temp_Cpl(iz)*thetazm(iz)
- ! endif
-
- ! if (temp_Fpl(iz)>0) then
- ! tmp = tmp + temp_Fpl(iz)*tempz(iz)
- ! else
- ! tmp = tmp + temp_Fpl(iz)*tempzm(iz)
- ! endif
-
- ! temp_rhs(iz,updatetlevel) = tmp &
- ! +temp_Hpl(iz)*vparzm(iz) & ! temp_Hpl is always negative
- ! +temp_Ipl(iz)*momentszm(3,iz) ! temp_Ipl is always negative
- ! end do
-
- case('fa2','fd4','we4')
- do iz = izs,ize
- temp_rhs(iz, updatetlevel) = diff_temp*tempzz(iz) & ! numerical diffusion added for stability
- +temp_Cpl(iz)*thetaz(iz) &
- +temp_Fpl(iz)*sqrt_exp_tempz(iz) &
- +temp_Hpl(iz)*vparz_n(iz) &
- +temp_Ipl(iz)*momentsz(3,iz)
-! - SQRT2*2._dp*nu*temp(iz,updatetlevel) ! Lenard-Bernstein Operator conserves energy
- end do
-
- end select
-
- ! Compute interpolation on v-grid, used by vpar_eq_rhs
- CALL interp_n2v(temp_rhs(:, updatetlevel), temp_rhs_v)
-
-END SUBROUTINE temp_eq_rhs
diff --git a/src/theta_eq_rhs.F90 b/src/theta_eq_rhs.F90
deleted file mode 100644
index 14ee49639b3d3a147489f4ee51cf72c759f669c7..0000000000000000000000000000000000000000
--- a/src/theta_eq_rhs.F90
+++ /dev/null
@@ -1,50 +0,0 @@
-SUBROUTINE theta_eq_rhs
-
- ! This can be used to evaluate the RHS of the Hermite Hierarchy Equation ....
-
- ! theta is ln(density)
- USE time_integration
- USE array
- USE fields
- USE space_grid
- USE model
-
-
- use prec_const
- IMPLICIT NONE
-
- INTEGER :: iz
- real(dp) :: tmp
-
-
- do iz = izs,ize ! Compute coefficients
- theta_Cpl(iz) = -sqrt_exp_temp(iz)*INVSQRT2*vpar_n(iz)
- theta_Hpl(iz) = -sqrt_exp_temp(iz)*SQRT2
- end do
-
-
- select case (gradient_scheme)
- ! case('up2') ! Upwind scheme, check the sign of coefficients
- ! do iz = izs,ize
- ! if (theta_Cpl(iz)>0) then
- ! tmp = theta_Cpl(iz)*thetaz(iz)
- ! else
- ! tmp = theta_Cpl(iz)*thetazm(iz)
- ! endif
-
- ! theta_rhs(iz, updatetlevel) = tmp &
- ! +theta_Hpl(iz)*vparzm(iz) ! theta_Hpl is always negative
-
- ! end do
-
- case('fa2','fd4','we4')
- do iz = izs,ize
- theta_rhs(iz, updatetlevel) = diff_theta*thetazz(iz) & ! numerical diffusion added for stability
- +theta_Cpl(iz)*thetaz(iz) &
- +theta_Hpl(iz)*vparz_n(iz)
-!theta_rhs(iz, updatetlevel) = 0._dp ! debug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- end do
-
- end select
-
-END SUBROUTINE theta_eq_rhs
diff --git a/src/utility_mod.F90 b/src/utility_mod.F90
index c22d1cc940118b61a7652968d790ded64a96b8bd..24445b854fb373b7f67c224c6e38eb67ce6fb7e2 100644
--- a/src/utility_mod.F90
+++ b/src/utility_mod.F90
@@ -53,18 +53,20 @@ CONTAINS
FUNCTION checkfield(field,str) RESULT(mlend)
- USE space_grid
+ USE fourier_grid
use prec_const
IMPLICIT NONE
- real(dp), DIMENSION(izs:ize), INTENT(IN) :: field
+ COMPLEX(dp), DIMENSION(ikrs:ikre,ikzs:ikze), INTENT(IN) :: field
CHARACTER(LEN=*), INTENT(IN) :: str
LOGICAL :: mlend
- real(dp) :: sumfield
+ COMPLEX(dp) :: sumfield
sumfield=SUM(field)
- mlend=is_nan(sumfield,str).OR.is_inf(sumfield,str)
+
+ mlend= is_nan( REAL(sumfield),str).OR.is_inf( REAL(sumfield),str) &
+ .OR. is_nan(AIMAG(sumfield),str).OR.is_inf(AIMAG(sumfield),str)
END FUNCTION checkfield
END MODULE utility
diff --git a/src/vpar_eq_rhs.F90 b/src/vpar_eq_rhs.F90
deleted file mode 100644
index c3b4aba9de914d9db809325dbaacafe64915fb5e..0000000000000000000000000000000000000000
--- a/src/vpar_eq_rhs.F90
+++ /dev/null
@@ -1,61 +0,0 @@
-SUBROUTINE vpar_eq_rhs
-
- USE time_integration
- USE array
- USE fields
- USE space_grid
- USE gradients
- USE model, ONLY: diff_vpar, gradient_scheme,nu
-
- use prec_const
- IMPLICIT NONE
-
- INTEGER :: iz
- real(dp) :: tmp
-
-
- do iz = izs,ize ! Compute coefficients
- vpar_Cpl(iz) = -sqrt_exp_temp_v(iz)*0.5_dp*INVSQRT2
- vpar_Dpl(iz) = INVSQRT2/sqrt_exp_temp_v(iz) ! ... should be +, otherwise instability
- vpar_Epl(iz) = -vpar(iz,updatetlevel)*0.5_dp
- vpar_Fpl(iz) = -(1._dp+2._dp*vpar(iz,updatetlevel)*vpar(iz,updatetlevel))*INVSQRT2
-! vpar_Fpl(iz) = -sqrt_exp_temp_v(iz)*(1._dp+2._dp*vpar(iz,updatetlevel)*vpar(iz,updatetlevel))*0.5_dp*INVSQRT2
- vpar_Hpl(iz) = -sqrt_exp_temp_v(iz)*SQRT2*vpar(iz,updatetlevel)
- end do
-
-
-
- select case (gradient_scheme)
- ! case('up2') ! Upwind scheme, check the sign of coefficients
- ! do iz = izs,ize
- ! if (vpar_Hpl(iz)>0) then
- ! tmp = vpar_Hpl(iz)*vparz(iz)
- ! else
- ! tmp = vpar_Hpl(iz)*vparzm(iz)
- ! endif
-
- ! vpar_rhs(iz,updatetlevel) = tmp &
- ! +vpar_Cpl(iz)*thetazm(iz) & ! vpar_Cpl is always negative
- ! +vpar_Dpl(iz)*phizm(iz) & ! vpar_Dpl is always negative
- ! +vpar_Epl(iz)*temp_rhs(iz,updatetlevel) &
- ! +vpar_Fpl(iz)*tempzm(iz) ! vpar_Fpl is always negative
- ! end do
-
- case('fa2','fd4','we4')
- do iz = izs,ize
- vpar_rhs(iz, updatetlevel) = diff_vpar*vparzz(iz) & ! numerical diffusion added for stability
- +vpar_Cpl(iz)*thetaz_v(iz) &
- +vpar_Dpl(iz)*phiz_v(iz) &
- +vpar_Epl(iz)*temp_rhs_v(iz) &
- +vpar_Fpl(iz)*sqrt_exp_tempz_v(iz) &
- +vpar_Hpl(iz)*vparz(iz) &
- -SQRT2*nu*vpar(iz,updatetlevel) ! add collisional friction from Lenard Bernstein Collision Operator
- end do
-
- end select
-
- ! Compute interpolation on n-grid, used by moments_eq_rhs
- CALL interp_v2n(vpar_rhs(:, updatetlevel), vpar_rhs_n)
-
-
-END SUBROUTINE vpar_eq_rhs
diff --git a/wk/fort.90 b/wk/fort.90
index 255f7ba9397e8d920c1ef4f90245fc58284b8746..966833a7245c764b8f7222f3bf1dc669ce9b526b 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -8,19 +8,22 @@ S. Brunner, T.M. Tran CRPP/EPFL
tmax=100 ! time normalized to 1/omega_pe
/
&GRID
- pmax=20 ! number of Hermite moments (including Ne,Te,Ue)
- nz=128
- zmin=0.
- zmax= 20 ! Normalized to lambda_De
+ pmaxe=10 ! number of Hermite moments (including Ne,Te,Ue)
+ jmaxe=5 ! number of Hermite moments (including Ne,Te,Ue)
+ pmaxi=10 ! number of Hermite moments (including Ne,Te,Ue)
+ jmaxi=5 ! number of Hermite moments (including Ne,Te,Ue)
+ nkr=1
+ krmin=0.
+ krmax=0. ! Normalized to ?
+ nkz=1
+ kzmin=1.
+ kzmax=1. ! Normalized to ?
/
&OUTPUT_PAR
nsave_0d=1000
- nsave_1d=1000
+ nsave_1d=0
nsave_2d=1000
nsave_3d=0
- write_theta=.true.
- write_temp=.true.
- write_vpar=.true.
write_moments=.true.
write_phi=.true.
write_doubleprecision=.true.
@@ -28,38 +31,24 @@ S. Brunner, T.M. Tran CRPP/EPFL
/
&MODEL_PAR
! Collisionality
- nu=0.01 ! Normalized collision frequency normalized to plasma frequency
- ! Numerically added diffusion
- diff_theta=0.001
- diff_temp=0.001
- diff_vpar=0.001
- diff_moments=0.001
- ! Other
- gradient_scheme='fd4'
- freeze_theta=.false.
- freeze_temp=.false.
- freeze_vpar=.false.
- freeze_moments=.false.
- freeze_phi=.false.
- fft_suppress_high_freq = .false.
- fft_sigma=1
+ nu = 0.01 ! Normalized collision frequency normalized to plasma frequency
+ tau_e = 1.0 ! T_e/T_e
+ tau_i = 1.0 ! T_i/T_e temperature ratio
+ sigma_e = 0.0233380 ! sqrt(m_e/m_i) mass ratio
+ sigma_i = 1.0 ! sqrt(m_i/m_i)
+ q_e =-1.0 ! Electrons charge
+ q_i = 1.0 ! Ions charge
+ eta_n = 1.0 ! L_perp / L_n Density gradient
+ eta_T = 0.0 ! L_perp / L_T Temperature gradient
+ eta_B = 0.5 ! L_perp / L_B Magnetic gradient and curvature
+ lambdaD = 0.0 ! Debye length
/
&INITIAL_CON
! Background value
- initback_density=1.
- initback_temp=1.
- initback_vpar=0.
initback_moments=0.
! Noise amplitude
- initnoise_density=0.
- initnoise_temp=0.
- initnoise_vpar=0.
initnoise_moments=0.
iseed=42
- init_nb_oscil_density=1. ! per domain, i.e., zmax
- init_nb_oscil_temp=1.
- init_ampli_density=1.0D-2
- init_ampli_temp=0
/
&TIME_INTEGRATION_PAR
numerical_scheme='RK4'