From b68fdf207edd9df676a40f2b08500a483bc9aed0 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Tue, 28 Mar 2023 17:07:43 +0200
Subject: [PATCH] Remove dp into xp which is meant to be set by user soon
---
src/array_mod.F90 | 62 ++++-----
src/basic_mod.F90 | 118 +++++++++---------
src/calculus_mod.F90 | 108 ++++++++--------
src/closure_mod.F90 | 16 +--
src/coeff_mod.F90 | 8 +-
src/collision_mod.F90 | 100 +++++++--------
src/control.F90 | 4 +-
src/cosolver_interface_mod.F90 | 38 +++---
src/diagnose.F90 | 38 +++---
src/fields_mod.F90 | 6 +-
src/fourier_mod.F90 | 6 +-
src/geometry_mod.F90 | 156 +++++++++++------------
src/ghosts_mod.F90 | 18 +--
src/grid_mod.F90 | 152 +++++++++++-----------
src/inital.F90 | 64 +++++-----
src/initial_par_mod.F90 | 6 +-
src/lag_interp_mod.F90 | 58 ++++-----
src/miller_mod.F90 | 136 ++++++++++----------
src/model_mod.F90 | 24 ++--
src/moments_eq_rhs_mod.F90 | 42 +++----
src/nonlinear_mod.F90 | 24 ++--
src/numerics_mod.F90 | 140 ++++++++++-----------
src/parallel_mod.F90 | 52 ++++----
src/prec_const_mod.F90 | 98 ++++++++-------
src/processing_mod.F90 | 136 ++++++++++----------
src/restarts_mod.F90 | 8 +-
src/solve_EM_fields.F90 | 20 +--
src/species_mod.F90 | 66 +++++-----
src/stepon.F90 | 8 +-
src/time_integration_mod.F90 | 222 ++++++++++++++++-----------------
src/utility_mod.F90 | 16 +--
wk/fast_analysis.m | 11 +-
32 files changed, 989 insertions(+), 972 deletions(-)
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index f4878080..74916783 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -4,65 +4,65 @@ MODULE array
implicit none
! Arrays to store the rhs, for time integration (ia,ip,ij,iky,ikx,iz)
- COMPLEX(dp), DIMENSION(:,:,:,:,:,:,:), ALLOCATABLE :: moments_rhs
+ COMPLEX(xp), DIMENSION(:,:,:,:,:,:,:), ALLOCATABLE :: moments_rhs
! Arrays of non-adiabatique moments
- COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: nadiab_moments
+ COMPLEX(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: nadiab_moments
! Derivatives and interpolated moments
- COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: ddz_napj
- COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: interp_napj
- COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: ddzND_Napj
+ COMPLEX(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: ddz_napj
+ COMPLEX(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: interp_napj
+ COMPLEX(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: ddzND_Napj
! Non linear term array (ip,ij,iky,ikx,iz)
- COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: Sapj ! electron
+ COMPLEX(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: Sapj ! electron
! a-a collision matrix (ia,ip,ij,iky,ikx,iz)
- REAL(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: Caa
+ REAL(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: Caa
! Test and field collision matrices (ia,ib,ip,ij,iky,ikx,iz)
- REAL(dp), DIMENSION(:,:,:,:,:,:,:), ALLOCATABLE :: Cab_F, Cab_T
+ REAL(xp), DIMENSION(:,:,:,:,:,:,:), ALLOCATABLE :: Cab_F, Cab_T
! nu x self collision matrix nuCself = nuaa*Caa + sum_b_neq_a nu_ab Cab_T
- REAL(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: nuCself
+ REAL(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: nuCself
! Collision term (ip,ij,iky,ikx,iz)
- COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: Capj
- COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: TColl_e_local, TColl_i_local
+ COMPLEX(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: Capj
+ COMPLEX(xp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: TColl_e_local, TColl_i_local
! dnjs coefficient storage (in, ij, is)
- COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: dnjs
+ COMPLEX(xp), DIMENSION(:,:,:), ALLOCATABLE :: dnjs
! Hermite fourth derivative coeff storage 4*sqrt(p!/(p-4)!)
- COMPLEX(dp), DIMENSION(:), ALLOCATABLE :: dv4_Hp_coeff
+ COMPLEX(xp), DIMENSION(:), ALLOCATABLE :: dv4_Hp_coeff
! lin rhs p,j coefficient storage (ip,ij)
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: xnapj
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xnapp1j, xnapp2j, xnapm1j, xnapm2j, xnapjp1, xnapjm1
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: ynapp1j, ynapm1j, ynapp1jm1, ynapm1jm1 ! mirror lin coeff for non adiab mom
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: zNapm1j, zNapm1jp1, zNapm1jm1 ! mirror lin coeff for adiab mom
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: xphij, xphijp1, xphijm1
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: xpsij, xpsijp1, xpsijm1
+ REAL(xp), DIMENSION(:,:,:), ALLOCATABLE :: xnapj
+ REAL(xp), DIMENSION(:,:), ALLOCATABLE :: xnapp1j, xnapp2j, xnapm1j, xnapm2j, xnapjp1, xnapjm1
+ REAL(xp), DIMENSION(:,:,:), ALLOCATABLE :: ynapp1j, ynapm1j, ynapp1jm1, ynapm1jm1 ! mirror lin coeff for non adiab mom
+ REAL(xp), DIMENSION(:,:,:), ALLOCATABLE :: zNapm1j, zNapm1jp1, zNapm1jm1 ! mirror lin coeff for adiab mom
+ REAL(xp), DIMENSION(:,:,:), ALLOCATABLE :: xphij, xphijp1, xphijm1
+ REAL(xp), DIMENSION(:,:,:), ALLOCATABLE :: xpsij, xpsijp1, xpsijm1
! Kernel function evaluation (ia,ij,iky,ikx,iz,odd/even p)
- REAL(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: kernel
+ REAL(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: kernel
! Poisson operator (iky,ikx,iz)
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_poisson_op
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_ampere_op
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_pol_ion
- ! REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: HF_phi_correction_operator
+ REAL(xp), DIMENSION(:,:,:), ALLOCATABLE :: inv_poisson_op
+ REAL(xp), DIMENSION(:,:,:), ALLOCATABLE :: inv_ampere_op
+ REAL(xp), DIMENSION(:,:,:), ALLOCATABLE :: inv_pol_ion
+ ! REAL(xp), DIMENSION(:,:,:), ALLOCATABLE :: HF_phi_correction_operator
! Kinetic spectrum sum_kx,ky(|Napj(z)|^2), (ia,ip,ij,iz) (should be real)
- REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Napjz
+ REAL(xp), DIMENSION(:,:,:,:), ALLOCATABLE :: Napjz
! particle density for electron and ions (iky,ikx,iz)
- COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: dens
+ COMPLEX(xp), DIMENSION(:,:,:,:), ALLOCATABLE :: dens
! particle fluid velocity for electron and ions (iky,ikx,iz)
- COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: upar
- COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: uper
+ COMPLEX(xp), DIMENSION(:,:,:,:), ALLOCATABLE :: upar
+ COMPLEX(xp), DIMENSION(:,:,:,:), ALLOCATABLE :: uper
! particle temperature for electron and ions (iky,ikx,iz)
- COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Tpar
- COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Tper
- COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: temp
+ COMPLEX(xp), DIMENSION(:,:,:,:), ALLOCATABLE :: Tpar
+ COMPLEX(xp), DIMENSION(:,:,:,:), ALLOCATABLE :: Tper
+ COMPLEX(xp), DIMENSION(:,:,:,:), ALLOCATABLE :: temp
END MODULE array
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index cc1c14f0..0784de19 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -1,18 +1,18 @@
MODULE basic
! Basic module for time dependent problems
use, intrinsic :: iso_c_binding
- use prec_const, ONLY : dp
+ use prec_const, ONLY : xp
IMPLICIT none
PRIVATE
! INCLUDE 'fftw3-mpi.f03'
! INPUT PARAMETERS
INTEGER, PUBLIC, PROTECTED :: nrun = 1 ! Number of time steps to run
- real(dp), PUBLIC, PROTECTED :: tmax = 100000.0 ! Maximum simulation time
- real(dp), PUBLIC, PROTECTED :: dt = 1.0 ! Time step
- real(dp), PUBLIC, PROTECTED :: maxruntime = 1e9 ! Maximum simulation CPU time
+ real(xp), PUBLIC, PROTECTED :: tmax = 100000.0 ! Maximum simulation time
+ real(xp), PUBLIC, PROTECTED :: dt = 1.0 ! Time step
+ real(xp), PUBLIC, PROTECTED :: maxruntime = 1e9 ! Maximum simulation CPU time
INTEGER, PUBLIC, PROTECTED :: job2load = 99 ! jobnum of the checkpoint to load
! Auxiliary variables
- real(dp), PUBLIC, PROTECTED :: time = 0 ! Current simulation time (Init from restart file)
+ real(xp), PUBLIC, PROTECTED :: time = 0 ! Current simulation time (Init from restart file)
INTEGER, PUBLIC, PROTECTED :: jobnum = 0 ! Job number
INTEGER, PUBLIC, PROTECTED :: step = 0 ! Calculation step of this run
@@ -32,9 +32,9 @@ MODULE basic
! To measure computation time
type :: chrono
- real(dp) :: tstart !start of the chrono
- real(dp) :: tstop !stop
- real(dp) :: ttot !cumulative time
+ real(xp) :: tstart !start of the chrono
+ real(xp) :: tstop !stop
+ real(xp) :: ttot !cumulative time
end type chrono
type(chrono), PUBLIC, PROTECTED :: chrono_runt, chrono_mrhs, chrono_advf, chrono_pois, chrono_sapj,&
@@ -47,8 +47,8 @@ MODULE basic
set_basic_cp, daytim, start_chrono, stop_chrono
INTERFACE allocate_array
- MODULE PROCEDURE allocate_array_dp1,allocate_array_dp2,allocate_array_dp3, &
- allocate_array_dp4, allocate_array_dp5, allocate_array_dp6, allocate_array_dp7
+ MODULE PROCEDURE allocate_array_xp1,allocate_array_xp2,allocate_array_xp3, &
+ allocate_array_xp4, allocate_array_xp5, allocate_array_xp6, allocate_array_xp7
MODULE PROCEDURE allocate_array_dc1,allocate_array_dc2,allocate_array_dc3, &
allocate_array_dc4, allocate_array_dc5, allocate_array_dc6, allocate_array_dc7
MODULE PROCEDURE allocate_array_i1,allocate_array_i2,allocate_array_i3,allocate_array_i4
@@ -56,7 +56,7 @@ MODULE basic
END INTERFACE allocate_array
INTERFACE str
- MODULE PROCEDURE str_dp, str_int
+ MODULE PROCEDURE str_xp, str_int
END INTERFACE
CONTAINS
@@ -125,7 +125,7 @@ CONTAINS
END SUBROUTINE
SUBROUTINE set_basic_cp(cstep_cp,time_cp,jobnum_cp)
IMPLICIT NONE
- REAL(dp), INTENT(IN) :: time_cp
+ REAL(xp), INTENT(IN) :: time_cp
INTEGER, INTENT(IN) :: cstep_cp, jobnum_cp
cstep = cstep_cp
time = time_cp
@@ -196,7 +196,7 @@ CONTAINS
SUBROUTINE display_h_min_s(time)
USE parallel, ONLY: my_id
IMPLICIT NONE
- real(dp) :: time
+ real(xp) :: time
integer :: days, hours, mins, secs
days = FLOOR(time/24./3600.);
hours= FLOOR(time/3600.);
@@ -228,13 +228,13 @@ CONTAINS
END SUBROUTINE display_h_min_s
!================================================================================
- function str_dp(k) result( str_ )
+ function str_xp(k) result( str_ )
! "Convert an integer to string."
- REAL(dp), intent(in) :: k
+ REAL(xp), intent(in) :: k
character(len=10):: str_
write (str_, "(G10.2)") k
str_ = adjustl(str_)
- end function str_dp
+ end function str_xp
function str_int(k) result( str_ )
! "Convert an integer to string."
@@ -245,117 +245,117 @@ CONTAINS
end function str_int
! To allocate arrays of doubles, integers, etc. at run time
- SUBROUTINE allocate_array_dp1(a,is1,ie1)
+ SUBROUTINE allocate_array_xp1(a,is1,ie1)
IMPLICIT NONE
- real(dp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
+ real(xp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1
ALLOCATE(a(is1:ie1))
- a=0.0_dp
- END SUBROUTINE allocate_array_dp1
+ a=0.0_xp
+ END SUBROUTINE allocate_array_xp1
- SUBROUTINE allocate_array_dp2(a,is1,ie1,is2,ie2)
+ SUBROUTINE allocate_array_xp2(a,is1,ie1,is2,ie2)
IMPLICIT NONE
- real(dp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ real(xp), 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
+ a=0.0_xp
+ END SUBROUTINE allocate_array_xp2
- SUBROUTINE allocate_array_dp3(a,is1,ie1,is2,ie2,is3,ie3)
+ SUBROUTINE allocate_array_xp3(a,is1,ie1,is2,ie2,is3,ie3)
IMPLICIT NONE
- real(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ real(xp), 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
+ a=0.0_xp
+ END SUBROUTINE allocate_array_xp3
- SUBROUTINE allocate_array_dp4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
+ SUBROUTINE allocate_array_xp4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
IMPLICIT NONE
- real(dp), DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ real(xp), 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
+ a=0.0_xp
+ END SUBROUTINE allocate_array_xp4
- SUBROUTINE allocate_array_dp5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
+ SUBROUTINE allocate_array_xp5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
IMPLICIT NONE
- real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ real(xp), 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=0.0_dp
- END SUBROUTINE allocate_array_dp5
+ a=0.0_xp
+ END SUBROUTINE allocate_array_xp5
- SUBROUTINE allocate_array_dp6(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5,is6,ie6)
+ SUBROUTINE allocate_array_xp6(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5,is6,ie6)
IMPLICIT NONE
- real(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ real(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5,is6,ie6
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5,is6:ie6))
- a=0.0_dp
- END SUBROUTINE allocate_array_dp6
+ a=0.0_xp
+ END SUBROUTINE allocate_array_xp6
- SUBROUTINE allocate_array_dp7(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5,is6,ie6,is7,ie7)
+ SUBROUTINE allocate_array_xp7(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5,is6,ie6,is7,ie7)
IMPLICIT NONE
- REAL(dp), DIMENSION(:,:,:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ REAL(xp), DIMENSION(:,:,:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5,is6,ie6,is7,ie7
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5,is6:ie6,is7:ie7))
- a=0.0_dp
- END SUBROUTINE allocate_array_dp7
+ a=0.0_xp
+ END SUBROUTINE allocate_array_xp7
!========================================
SUBROUTINE allocate_array_dc1(a,is1,ie1)
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
+ COMPLEX(xp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1
ALLOCATE(a(is1:ie1))
- a=CMPLX(0.0_dp,0.0_dp)
+ a=CMPLX(0.0_xp,0.0_xp)
END SUBROUTINE allocate_array_dc1
SUBROUTINE allocate_array_dc2(a,is1,ie1,is2,ie2)
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ COMPLEX(xp), 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)
+ a=CMPLX(0.0_xp,0.0_xp)
END SUBROUTINE allocate_array_dc2
SUBROUTINE allocate_array_dc3(a,is1,ie1,is2,ie2,is3,ie3)
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ COMPLEX(xp), 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)
+ a=CMPLX(0.0_xp,0.0_xp)
END SUBROUTINE allocate_array_dc3
SUBROUTINE allocate_array_dc4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ COMPLEX(xp), 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)
+ a=CMPLX(0.0_xp,0.0_xp)
END SUBROUTINE allocate_array_dc4
SUBROUTINE allocate_array_dc5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ COMPLEX(xp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5))
- a=CMPLX(0.0_dp,0.0_dp)
+ a=CMPLX(0.0_xp,0.0_xp)
END SUBROUTINE allocate_array_dc5
SUBROUTINE allocate_array_dc6(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5,is6,ie6)
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ COMPLEX(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5,is6,ie6
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5,is6:ie6))
- a=CMPLX(0.0_dp,0.0_dp)
+ a=CMPLX(0.0_xp,0.0_xp)
END SUBROUTINE allocate_array_dc6
SUBROUTINE allocate_array_dc7(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5,is6,ie6,is7,ie7)
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(:,:,:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ COMPLEX(xp), DIMENSION(:,:,:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5,is6,ie6,is7,ie7
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5,is6:ie6,is7:ie7))
- a=CMPLX(0.0_dp,0.0_dp)
+ a=CMPLX(0.0_xp,0.0_xp)
END SUBROUTINE allocate_array_dc7
!========================================
@@ -393,7 +393,7 @@ CONTAINS
SUBROUTINE allocate_array_i5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
IMPLICIT NONE
- real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
+ INTEGER, 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=0
diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
index a6c75d03..1e186bcb 100644
--- a/src/calculus_mod.F90
+++ b/src/calculus_mod.F90
@@ -1,21 +1,21 @@
MODULE calculus
! Routine to evaluate gradients, interpolation schemes and integrals
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
IMPLICIT NONE
- REAL(dp), dimension(-2:2) :: dz_usu = &
- (/ 1._dp/12._dp, -2._dp/3._dp, 0._dp, 2._dp/3._dp, -1._dp/12._dp /) ! fd4 centered stencil
- REAL(dp), dimension(-2:1) :: dz_o2e = &
- (/ 1._dp/24._dp,-9._dp/8._dp, 9._dp/8._dp,-1._dp/24._dp /) ! fd4 odd to even stencil
- REAL(dp), dimension(-1:2) :: dz_e2o = &
- (/ 1._dp/24._dp,-9._dp/8._dp, 9._dp/8._dp,-1._dp/24._dp /) ! fd4 odd to even stencil
- REAL(dp), dimension(-2:2) :: dz2_usu = &
- (/-1._dp/12._dp, 4._dp/3._dp, -5._dp/2._dp, 4._dp/3._dp, -1._dp/12._dp /)! 2th derivative, 4th order (for parallel hypdiff)
- REAL(dp), dimension(-2:2) :: dz4_usu = &
- (/ 1._dp, -4._dp, 6._dp, -4._dp, 1._dp /) ! 4th derivative, 2nd order (for parallel hypdiff)
- REAL(dp), dimension(-2:1) :: iz_o2e = &
- (/ -1._dp/16._dp, 9._dp/16._dp, 9._dp/16._dp, -1._dp/16._dp /) ! grid interpolation, 4th order, odd to even
- REAL(dp), dimension(-1:2) :: iz_e2o = &
- (/ -1._dp/16._dp, 9._dp/16._dp, 9._dp/16._dp, -1._dp/16._dp /) ! grid interpolation, 4th order, even to odd
+ REAL(xp), dimension(-2:2) :: dz_usu = &
+ (/ 1._xp/12._xp, -2._xp/3._xp, 0._xp, 2._xp/3._xp, -1._xp/12._xp /) ! fd4 centered stencil
+ REAL(xp), dimension(-2:1) :: dz_o2e = &
+ (/ 1._xp/24._xp,-9._xp/8._xp, 9._xp/8._xp,-1._xp/24._xp /) ! fd4 odd to even stencil
+ REAL(xp), dimension(-1:2) :: dz_e2o = &
+ (/ 1._xp/24._xp,-9._xp/8._xp, 9._xp/8._xp,-1._xp/24._xp /) ! fd4 odd to even stencil
+ REAL(xp), dimension(-2:2) :: dz2_usu = &
+ (/-1._xp/12._xp, 4._xp/3._xp, -5._xp/2._xp, 4._xp/3._xp, -1._xp/12._xp /)! 2th derivative, 4th order (for parallel hypdiff)
+ REAL(xp), dimension(-2:2) :: dz4_usu = &
+ (/ 1._xp, -4._xp, 6._xp, -4._xp, 1._xp /) ! 4th derivative, 2nd order (for parallel hypdiff)
+ REAL(xp), dimension(-2:1) :: iz_o2e = &
+ (/ -1._xp/16._xp, 9._xp/16._xp, 9._xp/16._xp, -1._xp/16._xp /) ! grid interpolation, 4th order, odd to even
+ REAL(xp), dimension(-1:2) :: iz_e2o = &
+ (/ -1._xp/16._xp, 9._xp/16._xp, 9._xp/16._xp, -1._xp/16._xp /) ! grid interpolation, 4th order, even to odd
PUBLIC :: simpson_rule_z, interp_z, grad_z, grad_z4, grad_z_5D, grad_z4_5D
CONTAINS
@@ -27,9 +27,9 @@ SUBROUTINE grad_z(target,local_nz,ngz,inv_deltaz,f,ddzf)
! not staggered : the derivative results must be on the same grid as the field
! staggered : the derivative is computed from a grid to the other
INTEGER, INTENT(IN) :: target, local_nz, ngz
- REAL(dp), INTENT(IN) :: inv_deltaz
- COMPLEX(dp),dimension(local_nz+ngz), INTENT(IN) :: f
- COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddzf
+ REAL(xp), INTENT(IN) :: inv_deltaz
+ COMPLEX(xp),dimension(local_nz+ngz), INTENT(IN) :: f
+ COMPLEX(xp),dimension(local_nz), INTENT(OUT) :: ddzf
INTEGER :: iz
IF(ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
SELECT CASE(TARGET)
@@ -46,16 +46,16 @@ SUBROUTINE grad_z(target,local_nz,ngz,inv_deltaz,f,ddzf)
ddzf(:) = ddzf(:) * inv_deltaz
END SELECT
ELSE
- ddzf = 0._dp
+ ddzf = 0._xp
ENDIF
CONTAINS
SUBROUTINE grad_z_o2e(local_nz,ngz,inv_deltaz,fo,ddzfe) ! Paruta 2018 eq (27)
! gives the gradient of a field from the odd grid to the even one
implicit none
INTEGER, INTENT(IN) :: local_nz, ngz
- REAL(dp), INTENT(IN) :: inv_deltaz
- COMPLEX(dp),dimension(local_nz+ngz), INTENT(IN) :: fo
- COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddzfe !
+ REAL(xp), INTENT(IN) :: inv_deltaz
+ COMPLEX(xp),dimension(local_nz+ngz), INTENT(IN) :: fo
+ COMPLEX(xp),dimension(local_nz), INTENT(OUT) :: ddzfe !
INTEGER :: iz
DO iz = 1,local_nz
ddzfe(iz) = dz_o2e(-2)*fo(iz ) + dz_o2e(-1)*fo(iz+1) &
@@ -68,9 +68,9 @@ CONTAINS
! gives the gradient of a field from the even grid to the odd one
implicit none
INTEGER, INTENT(IN) :: local_nz, ngz
- REAL(dp), INTENT(IN) :: inv_deltaz
- COMPLEX(dp),dimension(local_nz+ngz), INTENT(IN) :: fe
- COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddzfo
+ REAL(xp), INTENT(IN) :: inv_deltaz
+ COMPLEX(xp),dimension(local_nz+ngz), INTENT(IN) :: fe
+ COMPLEX(xp),dimension(local_nz), INTENT(OUT) :: ddzfo
INTEGER :: iz
DO iz = 1,local_nz
ddzfo(iz) = dz_e2o(-1)*fe(iz+1) + dz_e2o(0)*fe(iz+2) &
@@ -87,9 +87,9 @@ SUBROUTINE grad_z_5D(local_nz,ngz,inv_deltaz,f,ddzf)
! not staggered : the derivative results must be on the same grid as the field
! staggered : the derivative is computed from a grid to the other
INTEGER, INTENT(IN) :: local_nz, ngz
- REAL(dp), INTENT(IN) :: inv_deltaz
- COMPLEX(dp),dimension(:,:,:,:,:,:), INTENT(IN) :: f
- COMPLEX(dp),dimension(:,:,:,:,:,:), INTENT(OUT) :: ddzf
+ REAL(xp), INTENT(IN) :: inv_deltaz
+ COMPLEX(xp),dimension(:,:,:,:,:,:), INTENT(IN) :: f
+ COMPLEX(xp),dimension(:,:,:,:,:,:), INTENT(OUT) :: ddzf
INTEGER :: iz
IF(ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
DO iz = 1,local_nz
@@ -99,7 +99,7 @@ SUBROUTINE grad_z_5D(local_nz,ngz,inv_deltaz,f,ddzf)
+dz_usu( 1)*f(:,:,:,:,:,iz+3) + dz_usu( 2)*f(:,:,:,:,:,iz+4))
ENDDO
ELSE
- ddzf = 0._dp
+ ddzf = 0._xp
ENDIF
END SUBROUTINE grad_z_5D
@@ -108,9 +108,9 @@ SUBROUTINE grad_z2(local_nz,ngz,inv_deltaz,f,ddz2f)
! Compute the second order fourth derivative for periodic boundary condition
implicit none
INTEGER, INTENT(IN) :: local_nz, ngz
- REAL(dp), INTENT(IN) :: inv_deltaz
- COMPLEX(dp),dimension(local_nz+ngz), INTENT(IN) :: f
- COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddz2f
+ REAL(xp), INTENT(IN) :: inv_deltaz
+ COMPLEX(xp),dimension(local_nz+ngz), INTENT(IN) :: f
+ COMPLEX(xp),dimension(local_nz), INTENT(OUT) :: ddz2f
INTEGER :: iz
IF(ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
DO iz = 1,local_nz
@@ -119,7 +119,7 @@ SUBROUTINE grad_z2(local_nz,ngz,inv_deltaz,f,ddz2f)
+dz2_usu( 1)*f(iz+3) + dz2_usu( 2)*f(iz+4)
ENDDO
ELSE
- ddz2f = 0._dp
+ ddz2f = 0._xp
ENDIF
ddz2f = ddz2f * inv_deltaz**2
END SUBROUTINE grad_z2
@@ -128,9 +128,9 @@ SUBROUTINE grad_z4_5D(local_nz,ngz,inv_deltaz,f,ddz4f)
! Compute the second order fourth derivative for periodic boundary condition
implicit none
INTEGER, INTENT(IN) :: local_nz, ngz
- REAL(dp), INTENT(IN) :: inv_deltaz
- COMPLEX(dp),dimension(:,:,:,:,:,:), INTENT(IN) :: f
- COMPLEX(dp),dimension(:,:,:,:,:,:), INTENT(OUT) :: ddz4f
+ REAL(xp), INTENT(IN) :: inv_deltaz
+ COMPLEX(xp),dimension(:,:,:,:,:,:), INTENT(IN) :: f
+ COMPLEX(xp),dimension(:,:,:,:,:,:), INTENT(OUT) :: ddz4f
INTEGER :: iz
IF(ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
DO iz = 1,local_nz
@@ -140,7 +140,7 @@ SUBROUTINE grad_z4_5D(local_nz,ngz,inv_deltaz,f,ddz4f)
+dz4_usu( 1)*f(:,:,:,:,:,iz+3) + dz4_usu( 2)*f(:,:,:,:,:,iz+4))
ENDDO
ELSE
- ddz4f = 0._dp
+ ddz4f = 0._xp
ENDIF
END SUBROUTINE grad_z4_5D
@@ -148,9 +148,9 @@ SUBROUTINE grad_z4(local_nz,ngz,inv_deltaz,f,ddz4f)
! Compute the second order fourth derivative for periodic boundary condition
implicit none
INTEGER, INTENT(IN) :: local_nz, ngz
- REAL(dp), INTENT(IN) :: inv_deltaz
- COMPLEX(dp),dimension(local_nz+ngz), INTENT(IN) :: f
- COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddz4f
+ REAL(xp), INTENT(IN) :: inv_deltaz
+ COMPLEX(xp),dimension(local_nz+ngz), INTENT(IN) :: f
+ COMPLEX(xp),dimension(local_nz), INTENT(OUT) :: ddz4f
INTEGER :: iz
IF(ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
DO iz = 1,local_nz
@@ -159,7 +159,7 @@ SUBROUTINE grad_z4(local_nz,ngz,inv_deltaz,f,ddz4f)
+dz4_usu( 1)*f(iz+3) + dz4_usu( 2)*f(iz+4)
ENDDO
ELSE
- ddz4f = 0._dp
+ ddz4f = 0._xp
ENDIF
ddz4f = ddz4f * inv_deltaz**4
END SUBROUTINE grad_z4
@@ -172,8 +172,8 @@ SUBROUTINE interp_z(target,local_nz,ngz,f_in,f_out)
implicit none
INTEGER, INTENT(IN) :: local_nz, ngz
INTEGER, intent(in) :: target ! target grid : 0 for even grid, 1 for odd
- COMPLEX(dp),dimension(local_nz+ngz), INTENT(IN) :: f_in
- COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: f_out
+ COMPLEX(xp),dimension(local_nz+ngz), INTENT(IN) :: f_in
+ COMPLEX(xp),dimension(local_nz), INTENT(OUT) :: f_out
SELECT CASE(TARGET)
CASE(1) ! output on even grid
CALL interp_o2e_z(local_nz,ngz,f_in,f_out)
@@ -187,8 +187,8 @@ CONTAINS
! gives the value of a field from the odd grid to the even one
implicit none
INTEGER, INTENT(IN) :: local_nz, ngz
- COMPLEX(dp),dimension(local_nz+ngz), INTENT(IN) :: fo
- COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: fe
+ COMPLEX(xp),dimension(local_nz+ngz), INTENT(IN) :: fo
+ COMPLEX(xp),dimension(local_nz), INTENT(OUT) :: fe
INTEGER :: iz
! 4th order interp
DO iz = 1,local_nz
@@ -201,8 +201,8 @@ CONTAINS
! gives the value of a field from the even grid to the odd one
implicit none
INTEGER, INTENT(IN) :: local_nz, ngz
- COMPLEX(dp),dimension(local_nz+ngz), INTENT(IN) :: fe
- COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: fo
+ COMPLEX(xp),dimension(local_nz+ngz), INTENT(IN) :: fe
+ COMPLEX(xp),dimension(local_nz), INTENT(OUT) :: fo
INTEGER :: iz
! 4th order interp
DO iz = 1,local_nz
@@ -216,25 +216,25 @@ END SUBROUTINE interp_z
SUBROUTINE simpson_rule_z(local_nz,zweights_SR,f,intf)
! integrate f(z) over z using the simpon's rule. Assume periodic boundary conditions (f(ize+1) = f(izs))
!from molix BJ Frei
- USE prec_const, ONLY: dp, onethird
+ USE prec_const, ONLY: xp, onethird
USE parallel, ONLY: num_procs_z, rank_z, comm_z, manual_0D_bcast
USE mpi
implicit none
INTEGER, INTENT(IN) :: local_nz
- REAL(dp), dimension(local_nz), intent(in) :: zweights_SR
- complex(dp),dimension(local_nz), intent(in) :: f
- COMPLEX(dp), intent(out) :: intf
- COMPLEX(dp) :: buffer, local_int
+ REAL(xp), dimension(local_nz), intent(in) :: zweights_SR
+ complex(xp),dimension(local_nz), intent(in) :: f
+ COMPLEX(xp), intent(out) :: intf
+ COMPLEX(xp) :: buffer, local_int
INTEGER :: root, i_, iz, ierr
! Buil local sum using the weights of composite Simpson's rule
- local_int = 0._dp
+ local_int = 0._xp
DO iz = 1,local_nz
local_int = local_int + zweights_SR(iz)*f(iz)
ENDDO
buffer = local_int
root = 0
!Gather manually among the rank_z=0 processes and perform the sum
- intf = 0._dp
+ intf = 0._xp
IF (num_procs_z .GT. 1) THEN
!! Everyone sends its local_sum to root = 0
IF (rank_z .NE. root) THEN
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index cd6e231a..9180feeb 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -8,7 +8,7 @@ CONTAINS
! Positive Oob indices are approximated with a model
SUBROUTINE apply_closure_model
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
USE model, ONLY: CLOS
USE grid, ONLY: local_nj,ngj, jarray,&
local_np,ngp, parray, dmax
@@ -28,7 +28,7 @@ SUBROUTINE apply_closure_model
j: DO ij = 1,local_nj+ngj
p: DO ip = 1,local_np+ngp
IF ( parray(ip)+2*jarray(ij) .GT. dmax) THEN
- moments(ia,ip,ij,:,:,:,updatetlevel) = 0._dp
+ moments(ia,ip,ij,:,:,:,updatetlevel) = 0._xp
ENDIF
ENDDO p
ENDDO j
@@ -40,7 +40,7 @@ END SUBROUTINE apply_closure_model
! Positive Oob indices are approximated with a model
SUBROUTINE ghosts_upper_truncation
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
USE grid, ONLY: local_np,ngp,local_pmax, pmax,&
local_nj,ngj,local_jmax, jmax
USE fields, ONLY: moments
@@ -51,20 +51,20 @@ SUBROUTINE ghosts_upper_truncation
! applies only for the processes that evolve the highest moment degree
IF(local_jmax .GE. Jmax) THEN
DO ig = 1,ngj/2
- moments(:,:,local_nj+ngj/2+ig,:,:,:,updatetlevel) = 0._dp
+ moments(:,:,local_nj+ngj/2+ig,:,:,:,updatetlevel) = 0._xp
ENDDO
ENDIF
! applies only for the process that has largest p index
IF(local_pmax .GE. Pmax) THEN
DO ig = 1,ngp/2
- moments(:,local_np+ngp/2+ig,:,:,:,:,updatetlevel) = 0._dp
+ moments(:,local_np+ngp/2+ig,:,:,:,:,updatetlevel) = 0._xp
ENDDO
ENDIF
END SUBROUTINE ghosts_upper_truncation
! Negative OoB indices are 0
SUBROUTINE ghosts_lower_truncation
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
USE grid, ONLY: ngp,ngj,local_pmin,local_jmin
USE fields, ONLY: moments
USE time_integration, ONLY: updatetlevel
@@ -73,13 +73,13 @@ SUBROUTINE ghosts_lower_truncation
! zero truncation, An=0 for n<0
IF(local_jmin .EQ. 0) THEN
DO ig = 1,ngj/2
- moments(:,:,ig,:,:,:,updatetlevel) = 0._dp
+ moments(:,:,ig,:,:,:,updatetlevel) = 0._xp
ENDDO
ENDIF
! applies only for the process that has lowest p index
IF(local_pmin .EQ. 0) THEN
DO ig = 1,ngp/2
- moments(:,ig,:,:,:,:,updatetlevel) = 0._dp
+ moments(:,ig,:,:,:,:,updatetlevel) = 0._xp
ENDDO
ENDIF
diff --git a/src/coeff_mod.F90 b/src/coeff_mod.F90
index cfd32243..3c4c9125 100644
--- a/src/coeff_mod.F90
+++ b/src/coeff_mod.F90
@@ -44,11 +44,11 @@ CONTAINS
XOUT = TO_FM('0.0')
!
DO n1=0,n
- AUXN =Lpjl(REAL(m,dp)-0.5_dp,REAL(n,dp),REAL(n1,dp))
+ AUXN =Lpjl(REAL(m,xp)-0.5_xp,REAL(n,xp),REAL(n1,xp))
DO k1=0,k
- AUXK =Lpjl(-0.5_dp,REAL(k,dp),REAL(k1,dp))
+ AUXK =Lpjl(-0.5_xp,REAL(k,xp),REAL(k1,xp))
DO s1=0,s
- AUXS = Lpjl(-0.5_dp,REAL(s,dp),REAL(s1,dp))
+ AUXS = Lpjl(-0.5_xp,REAL(s,xp),REAL(s1,xp))
XOUT = XOUT + FACTORIAL(TO_FM(n1 + k1 + s1 ))*AUXN*AUXK*AUXS
ENDDO
ENDDO
@@ -75,7 +75,7 @@ CONTAINS
IMPLICIT NONE
!
- REAL(dp), intent(in) :: p,j,l
+ REAL(xp), intent(in) :: p,j,l
TYPE(FM) :: XOUT
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index a223da4d..903e6ddb 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -1,13 +1,13 @@
module collision
! contains the Hermite-Laguerre collision operators without the use of COSOLVER
-USE prec_const, ONLY : dp
+USE prec_const, ONLY : xp
IMPLICIT NONE
PRIVATE
CHARACTER(len=32), PUBLIC, PROTECTED :: collision_model ! (Lenard-Bernstein: 'LB', Dougherty: 'DG', Sugama: 'SG', Lorentz: 'LR', Landau: 'LD')
LOGICAL, PUBLIC, PROTECTED :: GK_CO =.true. ! activates GK effects on CO
LOGICAL, PUBLIC, PROTECTED :: INTERSPECIES =.true. ! activates interpecies collision
CHARACTER(len=128), PUBLIC, PROTECTED :: mat_file ! COSOlver matrix file names
-REAL(dp), PUBLIC, PROTECTED :: collision_kcut = 100.0
+REAL(xp), PUBLIC, PROTECTED :: collision_kcut = 100.0
LOGICAL, PUBLIC, PROTECTED :: cosolver_coll ! check if cosolver matrices are used
PUBLIC :: init_collision
@@ -96,12 +96,12 @@ CONTAINS
CASE ('SG','LR','LD')
CALL compute_cosolver_coll(GK_CO)
CASE ('none','hypcoll','dvpar4')
- Capj = 0._dp
+ Capj = 0._xp
CASE DEFAULT
ERROR STOP '>> ERROR << collision operator not recognized!!'
END SELECT
ELSE
- Capj = 0._dp
+ Capj = 0._xp
ENDIF
END SUBROUTINE compute_Capj
@@ -116,8 +116,8 @@ CONTAINS
USE array, ONLY: Capj
USE fields, ONLY: moments
IMPLICIT NONE
- COMPLEX(dp) :: TColl_
- REAL(dp) :: j_dp, p_dp, ba_2, kp
+ COMPLEX(xp) :: TColl_
+ REAL(xp) :: j_xp, p_xp, ba_2, kp
INTEGER :: iz,ikx,iky,ij,ip,ia,eo
DO iz = izs,ize
DO ikx = ikxs, ikxe;
@@ -127,16 +127,16 @@ CONTAINS
DO ia = ias, iae
!** Auxiliary variables **
eo = MODULO(parray(ip),2)
- p_dp = REAL(parray(ip),dp)
- j_dp = REAL(jarray(ij),dp)
+ p_xp = REAL(parray(ip),xp)
+ j_xp = REAL(jarray(ij),xp)
kp = kparray(iky,ikx,iz,eo)
ba_2 = kp**2 * sigma2_tau_o2(ia) ! this is (ba/2)^2
!** Assembling collison operator **
! -nuee (p + 2j) Nepj
- TColl_ = -nu_ab(ia,ia) * (p_dp + 2._dp*j_dp)*moments(ia,ip,ij,iky,ikx,iz,updatetlevel)
+ TColl_ = -nu_ab(ia,ia) * (p_xp + 2._xp*j_xp)*moments(ia,ip,ij,iky,ikx,iz,updatetlevel)
IF(GK_CO) THEN
- TColl_ = TColl_ - nu_ab(ia,ia) *2._dp*ba_2*moments(ia,ip,ij,iky,ikx,iz,updatetlevel)
+ TColl_ = TColl_ - nu_ab(ia,ia) *2._xp*ba_2*moments(ia,ip,ij,iky,ikx,iz,updatetlevel)
ENDIF
Capj(ia,ip,ij,iky,ikx,iz) = TColl_
ENDDO
@@ -170,11 +170,11 @@ CONTAINS
USE time_integration, ONLY: updatetlevel
USE array, ONLY: Capj
USE fields, ONLY: moments
- USE prec_const, ONLY: dp, SQRT2, twothird
+ USE prec_const, ONLY: xp, SQRT2, twothird
IMPLICIT NONE
- COMPLEX(dp) :: Tmp
+ COMPLEX(xp) :: Tmp
INTEGER :: iz,ikx,iky,ij,ip,ia, ipi,iji,izi
- REAL(dp) :: j_dp, p_dp
+ REAL(xp) :: j_xp, p_xp
DO iz = 1,local_nz
izi = iz + ngz/2
DO ikx = 1,local_nkx
@@ -185,17 +185,17 @@ CONTAINS
ipi = ip + ngp/2
DO ia = 1,local_na
!** Auxiliary variables **
- p_dp = REAL(parray(ipi),dp)
- j_dp = REAL(jarray(iji),dp)
+ p_xp = REAL(parray(ipi),xp)
+ j_xp = REAL(jarray(iji),xp)
!** Assembling collison operator **
- Tmp = -(p_dp + 2._dp*j_dp)*moments(ia,ipi,iji,iky,ikx,izi,updatetlevel)
- IF( (p_dp .EQ. 1._dp) .AND. (j_dp .EQ. 0._dp)) THEN !Ce10
+ Tmp = -(p_xp + 2._xp*j_xp)*moments(ia,ipi,iji,iky,ikx,izi,updatetlevel)
+ IF( (p_xp .EQ. 1._xp) .AND. (j_xp .EQ. 0._xp)) THEN !Ce10
Tmp = Tmp + moments(ia,ip1,ij1,iky,ikx,iz,updatetlevel)
- ELSEIF( (p_dp .EQ. 2._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! Ce20
+ ELSEIF( (p_xp .EQ. 2._xp) .AND. (j_xp .EQ. 0._xp)) THEN ! Ce20
Tmp = Tmp + twothird*moments(ia,ip2,ij0,iky,ikx,izi,updatetlevel) &
- SQRT2*twothird*moments(ia,ip0,ij1,iky,ikx,izi,updatetlevel)
- ELSEIF( (p_dp .EQ. 0._dp) .AND. (j_dp .EQ. 1._dp)) THEN ! Ce01
- Tmp = Tmp + 2._dp*twothird*moments(ia,ip0,ij1,iky,ikx,izi,updatetlevel) &
+ ELSEIF( (p_xp .EQ. 0._xp) .AND. (j_xp .EQ. 1._xp)) THEN ! Ce01
+ Tmp = Tmp + 2._xp*twothird*moments(ia,ip0,ij1,iky,ikx,izi,updatetlevel) &
-SQRT2*twothird*moments(ia,ip2,ij0,iky,ikx,izi,updatetlevel)
ENDIF
Capj(ia,ip,ij,iky,ikx,iz) = nu_ab(ia,ia) * Tmp
@@ -215,13 +215,13 @@ CONTAINS
local_nky, local_nkx, local_nz, ngz, kparray, nzgrid
USE species, ONLY: sigma2_tau_o2, sqrt_sigma2_tau_o2, nu_ab
USE array, ONLY: Capj, nadiab_moments, kernel
- USE prec_const, ONLY: dp, SQRT2, twothird
+ USE prec_const, ONLY: xp, SQRT2, twothird
IMPLICIT NONE
!! Local variables
- COMPLEX(dp) :: dens_,upar_,uperp_,Tpar_,Tperp_,Temp_
- COMPLEX(dp) :: nadiab_moment_0j, Tmp
- REAL(dp) :: Knp0, Knp1, Knm1, kp
- REAL(dp) :: n_dp, j_dp, p_dp, ba, ba_2
+ COMPLEX(xp) :: dens_,upar_,uperp_,Tpar_,Tperp_,Temp_
+ COMPLEX(xp) :: nadiab_moment_0j, Tmp
+ REAL(xp) :: Knp0, Knp1, Knm1, kp
+ REAL(xp) :: n_xp, j_xp, p_xp, ba, ba_2
INTEGER :: iz,ikx,iky,ij,ip,ia,eo,in, ipi,iji,izi,ini
DO iz = 1,local_nz
izi = iz + ngz/2
@@ -233,25 +233,25 @@ CONTAINS
ipi = ip + ngp/2
DO ia = 1,local_na
!** Auxiliary variables **
- p_dp = REAL(parray(ipi),dp)
- j_dp = REAL(jarray(iji),dp)
+ p_xp = REAL(parray(ipi),xp)
+ j_xp = REAL(jarray(iji),xp)
eo = MIN(nzgrid,MODULO(parray(ipi),2)+1)
kp = kparray(iky,ikx,izi,eo)
ba_2 = kp**2 * sigma2_tau_o2(ia) ! this is (l_a/2)^2
- ba = 2_dp*kp * sqrt_sigma2_tau_o2(ia) ! this is l_a
+ ba = 2_xp*kp * sqrt_sigma2_tau_o2(ia) ! this is l_a
!** Assembling collison operator **
! Velocity-space diffusion (similar to Lenard Bernstein)
! -nu (p + 2j + b^2/2) Napj
- Tmp = -(p_dp + 2._dp*j_dp + 2._dp*ba_2)*nadiab_moments(ia,ipi,iji,iky,ikx,izi)
+ Tmp = -(p_xp + 2._xp*j_xp + 2._xp*ba_2)*nadiab_moments(ia,ipi,iji,iky,ikx,izi)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
+ IF( p_xp .EQ. 0 ) THEN ! Kronecker p0
!** build required fluid moments **
- dens_ = 0._dp
- upar_ = 0._dp; uperp_ = 0._dp
- Tpar_ = 0._dp; Tperp_ = 0._dp
+ dens_ = 0._xp
+ upar_ = 0._xp; uperp_ = 0._xp
+ Tpar_ = 0._xp; Tperp_ = 0._xp
DO in = 1,local_nj
ini = in + ngj/2
- n_dp = REAL(jarray(ini),dp)
+ n_xp = REAL(jarray(ini),xp)
! Store the kernels for sparing readings
Knp0 = kernel(ia,ini ,iky,ikx,izi,eo)
Knp1 = kernel(ia,ini+1,iky,ikx,izi,eo)
@@ -261,40 +261,40 @@ CONTAINS
! Density
dens_ = dens_ + Knp0 * nadiab_moment_0j
! Perpendicular velocity
- uperp_ = uperp_ + ba*0.5_dp*(Knp0 - Knm1) * nadiab_moment_0j
+ uperp_ = uperp_ + ba*0.5_xp*(Knp0 - Knm1) * nadiab_moment_0j
! Parallel temperature
Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments(ia,ip2,ini,iky,ikx,izi) + nadiab_moment_0j)
! Perpendicular temperature
- Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
+ Tperp_ = Tperp_ + ((2._xp*n_xp+1._xp)*Knp0 - (n_xp+1._xp)*Knp1 - n_xp*Knm1)*nadiab_moment_0j
ENDDO
- Temp_ = (Tpar_ + 2._dp*Tperp_)/3._dp - dens_
+ Temp_ = (Tpar_ + 2._xp*Tperp_)/3._xp - dens_
! Add energy restoring term
- Tmp = Tmp + Temp_* 4._dp * j_dp * kernel(ia,iji ,iky,ikx,izi,eo)
- Tmp = Tmp - Temp_* 2._dp * (j_dp + 1._dp) * kernel(ia,iji+1,iky,ikx,izi,eo)
- Tmp = Tmp - Temp_* 2._dp * j_dp * kernel(ia,iji-1,iky,ikx,izi,eo)
+ Tmp = Tmp + Temp_* 4._xp * j_xp * kernel(ia,iji ,iky,ikx,izi,eo)
+ Tmp = Tmp - Temp_* 2._xp * (j_xp + 1._xp) * kernel(ia,iji+1,iky,ikx,izi,eo)
+ Tmp = Tmp - Temp_* 2._xp * j_xp * kernel(ia,iji-1,iky,ikx,izi,eo)
Tmp = Tmp + uperp_*ba* (kernel(ia,iji,iky,ikx,izi,eo) - kernel(ia,iji-1,iky,ikx,izi,eo))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ELSEIF( p_dp .eq. 1 ) THEN ! kronecker p1
+ ELSEIF( p_xp .eq. 1 ) THEN ! kronecker p1
!** build required fluid moments **
- upar_ = 0._dp
+ upar_ = 0._xp
DO in = 1,local_nj
ini = in + ngj/2
- n_dp = REAL(jarray(ini),dp)
+ n_xp = REAL(jarray(ini),xp)
! Parallel velocity
upar_ = upar_ + kernel(ia,ini,iky,ikx,izi,eo) * nadiab_moments(ia,ip1,ini,iky,ikx,izi)
ENDDO
Tmp = Tmp + upar_*kernel(ia,iji,iky,ikx,izi,eo)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ELSEIF( p_dp .eq. 2 ) THEN ! kronecker p2
+ ELSEIF( p_xp .eq. 2 ) THEN ! kronecker p2
!** build required fluid moments **
- dens_ = 0._dp
- upar_ = 0._dp; uperp_ = 0._dp
- Tpar_ = 0._dp; Tperp_ = 0._dp
+ dens_ = 0._xp
+ upar_ = 0._xp; uperp_ = 0._xp
+ Tpar_ = 0._xp; Tperp_ = 0._xp
DO in = 1,local_nj
ini = in + ngj/2
- n_dp = REAL(jarray(ini),dp)
+ n_xp = REAL(jarray(ini),xp)
! Store the kernels for sparing readings
Knp0 = kernel(ia,ini ,iky,ikx,izi,eo)
Knp1 = kernel(ia,ini+1,iky,ikx,izi,eo)
@@ -306,9 +306,9 @@ CONTAINS
! Parallel temperature
Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments(ia,ip2,ini,iky,ikx,izi) + nadiab_moment_0j)
! Perpendicular temperature
- Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
+ Tperp_ = Tperp_ + ((2._xp*n_xp+1._xp)*Knp0 - (n_xp+1._xp)*Knp1 - n_xp*Knm1)*nadiab_moment_0j
ENDDO
- Temp_ = (Tpar_ + 2._dp*Tperp_)/3._dp - dens_
+ Temp_ = (Tpar_ + 2._xp*Tperp_)/3._xp - dens_
Tmp = Tmp + Temp_*SQRT2*kernel(ia,iji,iky,ikx,izi,eo)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ENDIF
diff --git a/src/control.F90 b/src/control.F90
index b70961cd..c8e307da 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -4,11 +4,11 @@ SUBROUTINE control
use basic, ONLY: str,daytim,speak,basic_data,&
nlend,step,increase_step,increase_time,increase_cstep,&
chrono_runt,chrono_step, chrono_diag, start_chrono, stop_chrono
- use prec_const, ONLY: dp, stdout
+ use prec_const, ONLY: xp, stdout
USE parallel, ONLY: ppinit
USE mpi
IMPLICIT NONE
- REAL(dp) :: t_init_diag_0, t_init_diag_1
+ REAL(xp) :: t_init_diag_0, t_init_diag_1
INTEGER :: ierr
! start the chronometer for the total runtime
CALL start_chrono(chrono_runt)
diff --git a/src/cosolver_interface_mod.F90 b/src/cosolver_interface_mod.F90
index f5b765e9..ade4d052 100644
--- a/src/cosolver_interface_mod.F90
+++ b/src/cosolver_interface_mod.F90
@@ -1,6 +1,6 @@
module cosolver_interface
! contains the Hermite-Laguerre collision operators solved using COSOlver.
-USE prec_const, ONLY: dp
+USE prec_const, ONLY: xp
IMPLICIT NONE
PRIVATE
PUBLIC :: load_COSOlver_mat, compute_cosolver_coll
@@ -19,9 +19,9 @@ CONTAINS
USE MPI
IMPLICIT NONE
LOGICAL, INTENT(IN) :: GK_CO
- COMPLEX(dp), DIMENSION(total_np) :: local_coll, buffer
- COMPLEX(dp), DIMENSION(local_np) :: TColl_distr
- COMPLEX(dp) :: Tmp_
+ COMPLEX(xp), DIMENSION(total_np) :: local_coll, buffer
+ COMPLEX(xp), DIMENSION(local_np) :: TColl_distr
+ COMPLEX(xp) :: Tmp_
INTEGER :: iz,ikx,iky,ij,ip,ia,ikx_C,iky_C,iz_C
INTEGER :: ierr
z:DO iz = 1,local_nz
@@ -74,7 +74,7 @@ CONTAINS
USE species, ONLY: nu_ab
IMPLICIT NONE
INTEGER, INTENT(IN) :: ia,ip,ij,iky,ikx,iz,ikx_C,iky_C,iz_C
- COMPLEX(dp), INTENT(OUT) :: local_coll
+ COMPLEX(xp), INTENT(OUT) :: local_coll
INTEGER :: ib,iq,il
INTEGER :: p_int,q_int,j_int,l_int
INTEGER :: izi, iqi, ili
@@ -82,7 +82,7 @@ CONTAINS
p_int = parray_full(ip)
j_int = jarray_full(ij)
!! Apply the cosolver collision matrix
- local_coll = 0._dp ! Initialization
+ local_coll = 0._xp ! Initialization
q:DO iq = 1,local_np
iqi = iq + ngp/2
q_int = parray(iqi)
@@ -124,20 +124,20 @@ CONTAINS
! Input
LOGICAL, INTENT(IN) :: GK_CO, INTERSPECIES
CHARACTER(len=128), INTENT(IN) :: matfile ! COSOlver matrix file names
- REAL(dp), INTENT(IN) :: collision_kcut
+ REAL(xp), INTENT(IN) :: collision_kcut
! Local variables
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Caa_full,CabT_full, CabF_full ! To load the self entire matrices
- REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Caa__kp ! To store the coeff that will be used along kperp
- REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: CabF_kp,CabT_kp ! ''
- REAL(dp), DIMENSION(:), ALLOCATABLE :: kp_grid_mat ! kperp grid of the matrices
- REAL(dp), DIMENSION(2) :: dims
+ REAL(xp), DIMENSION(:,:), ALLOCATABLE :: Caa_full,CabT_full, CabF_full ! To load the self entire matrices
+ REAL(xp), DIMENSION(:,:,:,:), ALLOCATABLE :: Caa__kp ! To store the coeff that will be used along kperp
+ REAL(xp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: CabF_kp,CabT_kp ! ''
+ REAL(xp), DIMENSION(:), ALLOCATABLE :: kp_grid_mat ! kperp grid of the matrices
+ REAL(xp), DIMENSION(2) :: dims
! Indices for row and columns of the COSOlver matrix (4D compressed 2D matrices)
INTEGER :: irow_sub, irow_full, icol_sub, icol_full
INTEGER :: fid ! file indexation
INTEGER :: ip, ij, il, ikx, iky, iz, ik, ikp, ikps_C,ikpe_C,ia,ib ! indices for loops
INTEGER :: pdim, jdim ! dimensions of the COSOlver matrices
INTEGER :: ikp_next, ikp_prev, nkp_mat, ikp_mat
- REAL(dp) :: kp_max,kperp_sim, kperp_mat, zerotoone
+ REAL(xp) :: kp_max,kperp_sim, kperp_mat, zerotoone
CHARACTER(len=128) :: var_name, ikp_string, name_a, name_b
CHARACTER(len=1) :: letter_a, letter_b
! Opening the compiled cosolver matrices results
@@ -225,7 +225,7 @@ CONTAINS
ENDDO
ENDDO
ELSE
- CabT_kp(ia,ib,:,:,:) = 0._dp; CabF_kp(ia,ib,:,:,:) = 0._dp
+ CabT_kp(ia,ib,:,:,:) = 0._xp; CabF_kp(ia,ib,:,:,:) = 0._xp
ENDIF
ENDDO b
ENDDO a
@@ -259,15 +259,15 @@ CONTAINS
! write(*,*) kp_grid_mat(ikp_prev), '<', kperp_sim, '<', kp_grid_mat(ikp_next)
ENDIF
! 0->1 variable for linear interp, i.e. zero2one = (k-k0)/(k1-k0)
- zerotoone = MIN(1._dp,(kperp_sim - kp_grid_mat(ikp_prev))/(kp_grid_mat(ikp_next) - kp_grid_mat(ikp_prev)))
+ zerotoone = MIN(1._xp,(kperp_sim - kp_grid_mat(ikp_prev))/(kp_grid_mat(ikp_next) - kp_grid_mat(ikp_prev)))
! Linear interpolation between previous and next kperp matrix values
Caa (:,:,:,iky,ikx,iz) = (Caa__kp(:,:,:,ikp_next) - Caa__kp(:,:,:,ikp_prev))*zerotoone + Caa__kp(:,:,:,ikp_prev)
IF(INTERSPECIES) THEN
Cab_T(:,:,:,:,iky,ikx,iz) = (CabT_kp(:,:,:,:,ikp_next) - CabT_kp(:,:,:,:,ikp_prev))*zerotoone + CabT_kp(:,:,:,:,ikp_prev)
Cab_F(:,:,:,:,iky,ikx,iz) = (CabF_kp(:,:,:,:,ikp_next) - CabF_kp(:,:,:,:,ikp_prev))*zerotoone + CabF_kp(:,:,:,:,ikp_prev)
ELSE
- Cab_T(:,:,:,:,iky,ikx,iz) = 0._dp
- Cab_F(:,:,:,:,iky,ikx,iz) = 0._dp
+ Cab_T(:,:,:,:,iky,ikx,iz) = 0._xp
+ Cab_F(:,:,:,:,iky,ikx,iz) = 0._xp
ENDIF
ENDDO
ENDDO
@@ -282,8 +282,8 @@ CONTAINS
IF( .NOT. INTERSPECIES ) THEN
CALL speak("--Like Species operator--")
- Cab_F = 0._dp;
- Cab_T = 0._dp;
+ Cab_F = 0._xp;
+ Cab_T = 0._xp;
ENDIF
! Build the self matrix
! nuCself = nuaa*Caa + sum_b_neq_a nu_ab Cab_T
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 78243c2a..3cf0cb9e 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -132,7 +132,7 @@ SUBROUTINE diagnose_full(kstep)
CALL putarr(fidres, "/data/grid/coordkx", kxarray_full, "kx*rho_s0", ionode=0)
CALL putarr(fidres, "/data/grid/coordky", kyarray_full, "ky*rho_s0", ionode=0)
CALL putarr(fidres, "/data/grid/coordz", zarray_full, "z/R", ionode=0)
- CALL putarr(fidres, "/data/grid/coordp" , parray_full, "p", ionode=0)
+ CALL putarr(fidres, "/data/grid/coorxp" , parray_full, "p", ionode=0)
CALL putarr(fidres, "/data/grid/coordj" , jarray_full, "j", ionode=0)
! Metric info
CALL creatg(fidres, "/data/metric", "Metric data")
@@ -289,7 +289,7 @@ SUBROUTINE diagnose_0d
CALL append(fidres, "/profiler/time", time,ionode=0)
! Processing data
CALL append(fidres, "/data/var0d/time", time,ionode=0)
- CALL append(fidres, "/data/var0d/cstep", real(cstep,dp),ionode=0)
+ CALL append(fidres, "/data/var0d/cstep", real(cstep,xp),ionode=0)
CALL getatt(fidres, "/data/var0d/", "frames",iframe0d)
iframe0d=iframe0d+1
CALL attach(fidres,"/data/var0d/" , "frames", iframe0d)
@@ -330,12 +330,12 @@ SUBROUTINE diagnose_3d
IMPLICIT NONE
CHARACTER :: letter_a
INTEGER :: ia
- COMPLEX(dp), DIMENSION(local_nky,local_nkx,local_nz) :: Na00_
- COMPLEX(dp), DIMENSION(local_nky,local_nkx,local_nz) :: fmom
- COMPLEX(dp), DIMENSION(local_np, local_nj, local_nz) :: Napjz_
+ COMPLEX(xp), DIMENSION(local_nky,local_nkx,local_nz) :: Na00_
+ COMPLEX(xp), DIMENSION(local_nky,local_nkx,local_nz) :: fmom
+ COMPLEX(xp), DIMENSION(local_np, local_nj, local_nz) :: Napjz_
! add current time, cstep and frame
CALL append(fidres, "/data/var3d/time", time,ionode=0)
- CALL append(fidres, "/data/var3d/cstep", real(cstep,dp),ionode=0)
+ CALL append(fidres, "/data/var3d/cstep", real(cstep,xp),ionode=0)
CALL getatt(fidres, "/data/var3d/", "frames",iframe3d)
iframe3d=iframe3d+1
CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
@@ -350,7 +350,7 @@ SUBROUTINE diagnose_3d
! gyrocenter density
Na00_ = moments(ia,ip0,ij0,:,:,(1+ngz/2):(local_nz+ngz/2),updatetlevel)
ELSE
- Na00_ = 0._dp
+ Na00_ = 0._xp
ENDIF
CALL write_field3d_kykxz(Na00_, 'N'//letter_a//'00')
! <<Napj>x>y spectrum
@@ -389,9 +389,9 @@ SUBROUTINE diagnose_3d
SUBROUTINE write_field3d_kykxz(field, text)
USE parallel, ONLY : gather_xyz, num_procs
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(local_nky,total_nkx,local_nz), INTENT(IN) :: field
+ COMPLEX(xp), DIMENSION(local_nky,total_nkx,local_nz), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
- COMPLEX(dp), DIMENSION(total_nky,total_nkx,total_nz) :: field_full
+ COMPLEX(xp), DIMENSION(total_nky,total_nkx,total_nz) :: field_full
CHARACTER(50) :: dset_name
field_full = 0;
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
@@ -407,9 +407,9 @@ SUBROUTINE diagnose_3d
SUBROUTINE write_field3d_pjz(field, text)
USE parallel, ONLY : gather_pjz, num_procs
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(local_np,local_nj,local_nz), INTENT(IN) :: field
+ COMPLEX(xp), DIMENSION(local_np,local_nj,local_nz), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
- COMPLEX(dp), DIMENSION(total_np,total_nj,total_nz) :: field_full
+ COMPLEX(xp), DIMENSION(total_np,total_nj,total_nz) :: field_full
CHARACTER(LEN=50) :: dset_name
field_full = 0;
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
@@ -434,11 +434,11 @@ SUBROUTINE diagnose_5d
ngp, ngj, ngz, total_na
USE time_integration, ONLY: updatetlevel, ntimelevel
USE diagnostics_par
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
IMPLICIT NONE
CALL append(fidres, "/data/var5d/time", time,ionode=0)
- CALL append(fidres, "/data/var5d/cstep", real(cstep,dp),ionode=0)
+ CALL append(fidres, "/data/var5d/cstep", real(cstep,xp),ionode=0)
CALL getatt(fidres, "/data/var5d/", "frames",iframe5d)
iframe5d=iframe5d+1
CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
@@ -453,12 +453,12 @@ SUBROUTINE diagnose_5d
USE basic, ONLY: GATHERV_OUTPUT, jobnum, dt
USE futils, ONLY: attach, putarr, putarrnd
USE parallel, ONLY: gather_pjxyz, num_procs
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(total_na,local_np+ngp,local_nj+ngj,local_nky,local_nkx,local_nz+ngz,ntimelevel), INTENT(IN) :: field
+ COMPLEX(xp), DIMENSION(total_na,local_np+ngp,local_nj+ngj,local_nky,local_nkx,local_nz+ngz,ntimelevel), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
- COMPLEX(dp), DIMENSION(total_na,local_np,local_nj,local_nky,local_nkx,local_nz) :: field_sub
- COMPLEX(dp), DIMENSION(total_na,total_np,total_nj,total_nky,total_nkx,total_nz) :: field_full
+ COMPLEX(xp), DIMENSION(total_na,local_np,local_nj,local_nky,local_nkx,local_nz) :: field_sub
+ COMPLEX(xp), DIMENSION(total_na,total_np,total_nj,total_nky,total_nkx,total_nz) :: field_full
CHARACTER(LEN=50) :: dset_name
field_sub = field(1:total_na,(1+ngp/2):(local_np+ngp/2),((1+ngj/2)):((local_nj+ngj/2)),&
1:local_nky,1:local_nkx, ((1+ngz/2)):((local_nz+ngz/2)),updatetlevel)
@@ -486,12 +486,12 @@ SUBROUTINE spit_snapshot_check
local_nky,local_nz, ngz
USE parallel, ONLY: gather_xyz, my_id
USE basic
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
IMPLICIT NONE
LOGICAL :: file_exist
INTEGER :: fid_check, ikx, iky, iz
CHARACTER(256) :: check_filename
- COMPLEX(dp), DIMENSION(total_nky,total_nkx,total_nz) :: field_to_check
+ COMPLEX(xp), DIMENSION(total_nky,total_nkx,total_nz) :: field_to_check
!! Spit a snapshot of PHI if requested (triggered by creating a file named "check_phi")
INQUIRE(file='check_phi', exist=file_exist)
IF( file_exist ) THEN
diff --git a/src/fields_mod.F90 b/src/fields_mod.F90
index ba3cc5cf..f689f8f5 100644
--- a/src/fields_mod.F90
+++ b/src/fields_mod.F90
@@ -4,13 +4,13 @@ MODULE fields
implicit none
!------------------MOMENTS Napj------------------
! Hermite-Moments: N_a^pj ! dimensions correspond to: species (a), p, j, kx, ky, z, updatetlevel.
- COMPLEX(dp), DIMENSION(:,:,:,:,:,:,:), ALLOCATABLE :: moments
+ COMPLEX(xp), DIMENSION(:,:,:,:,:,:,:), ALLOCATABLE :: moments
!------------------ELECTROSTATIC POTENTIAL------------------
! Normalized electric potential: \hat{\phi} ! (kx,ky,z)
- COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: phi
+ COMPLEX(xp), DIMENSION(:,:,:), ALLOCATABLE :: phi
!------------------Vector field part
- COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: psi
+ COMPLEX(xp), DIMENSION(:,:,:), ALLOCATABLE :: psi
END MODULE fields
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index f6e226c2..28acf05f 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -72,9 +72,9 @@ MODULE fourier
local_nky_ptr, local_nkx_ptr, F_, G_, sum_real_)
IMPLICIT NONE
INTEGER(C_INTPTR_T), INTENT(IN) :: local_nkx_ptr,local_nky_ptr
- REAL(dp), INTENT(IN) :: inv_Nx, inv_Ny
- REAL(dp), DIMENSION(local_nky_ptr), INTENT(IN) :: ky_, AA_y
- REAL(dp), DIMENSION(local_nkx_ptr), INTENT(IN) :: kx_, AA_x
+ REAL(xp), INTENT(IN) :: inv_Nx, inv_Ny
+ REAL(xp), DIMENSION(local_nky_ptr), INTENT(IN) :: ky_, AA_y
+ REAL(xp), DIMENSION(local_nkx_ptr), INTENT(IN) :: kx_, AA_x
COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(local_nky_ptr,local_nkx_ptr), &
INTENT(IN) :: F_(:,:), G_(:,:)
real(C_DOUBLE), pointer, INTENT(INOUT) :: sum_real_(:,:)
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index c354536f..1ee9153c 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -1,25 +1,25 @@
module geometry
! computes geometrical quantities
! Adapted from B.J.Frei MOLIX code (2021)
- use prec_const, ONLY: dp
+ use prec_const, ONLY: xp
implicit none
PRIVATE
! Geometry input parameters
CHARACTER(len=16), &
PUBLIC, PROTECTED :: geom = 's-alpha'
- REAL(dp), PUBLIC, PROTECTED :: q0 = 1.4_dp ! safety factor
- REAL(dp), PUBLIC, PROTECTED :: shear = 0._dp ! magnetic field shear
- REAL(dp), PUBLIC, PROTECTED :: eps = 0.18_dp ! inverse aspect ratio
- REAL(dp), PUBLIC, PROTECTED :: alpha_MHD = 0 ! shafranov shift effect alpha = -q2 R dbeta/dr
+ REAL(xp), PUBLIC, PROTECTED :: q0 = 1.4_xp ! safety factor
+ REAL(xp), PUBLIC, PROTECTED :: shear = 0._xp ! magnetic field shear
+ REAL(xp), PUBLIC, PROTECTED :: eps = 0.18_xp ! inverse aspect ratio
+ REAL(xp), PUBLIC, PROTECTED :: alpha_MHD = 0 ! shafranov shift effect alpha = -q2 R dbeta/dr
! parameters for Miller geometry
- REAL(dp), PUBLIC, PROTECTED :: kappa = 1._dp ! elongation (1 for circular)
- REAL(dp), PUBLIC, PROTECTED :: s_kappa = 0._dp ! r normalized derivative skappa = r/kappa dkappa/dr
- REAL(dp), PUBLIC, PROTECTED :: delta = 0._dp ! triangularity
- REAL(dp), PUBLIC, PROTECTED :: s_delta = 0._dp ! '' sdelta = r/sqrt(1-delta2) ddelta/dr
- REAL(dp), PUBLIC, PROTECTED :: zeta = 0._dp ! squareness
- REAL(dp), PUBLIC, PROTECTED :: s_zeta = 0._dp ! '' szeta = r dzeta/dr
+ REAL(xp), PUBLIC, PROTECTED :: kappa = 1._xp ! elongation (1 for circular)
+ REAL(xp), PUBLIC, PROTECTED :: s_kappa = 0._xp ! r normalized derivative skappa = r/kappa dkappa/dr
+ REAL(xp), PUBLIC, PROTECTED :: delta = 0._xp ! triangularity
+ REAL(xp), PUBLIC, PROTECTED :: s_delta = 0._xp ! '' sdelta = r/sqrt(1-delta2) ddelta/dr
+ REAL(xp), PUBLIC, PROTECTED :: zeta = 0._xp ! squareness
+ REAL(xp), PUBLIC, PROTECTED :: s_zeta = 0._xp ! '' szeta = r dzeta/dr
! to apply shift in the parallel z-BC if shearless
- REAL(dp), PUBLIC, PROTECTED :: shift_y = 0._dp ! for Arno <3
+ REAL(xp), PUBLIC, PROTECTED :: shift_y = 0._xp ! for Arno <3
INTEGER, PUBLIC, PROTECTED :: Npol = 1 ! number of poloidal turns
! Chooses the type of parallel BC we use for the unconnected kx modes (active for non-zero shear only)
! 'periodic' : Connect a disconnected kx to a mode on the other cadran
@@ -30,39 +30,39 @@ implicit none
PUBLIC, PROTECTED :: parallel_bc
! GENE unused additional parameters for miller_mod
- REAL(dp), PUBLIC, PROTECTED :: edge_opt = 0._dp ! meant to redistribute the points in z
- REAL(dp), PUBLIC, PROTECTED :: major_R = 1._dp ! major radius
- REAL(dp), PUBLIC, PROTECTED :: major_Z = 0._dp ! vertical elevation
- REAL(dp), PUBLIC, PROTECTED :: dpdx_pm_geom = 0._dp ! amplitude mag. eq. pressure grad.
- REAL(dp), PUBLIC, PROTECTED :: C_y = 0._dp ! defines y coordinate : Cy (q theta - phi)
- REAL(dp), PUBLIC, PROTECTED :: C_xy = 1._dp ! defines x coordinate : B = Cxy Vx x Vy
+ REAL(xp), PUBLIC, PROTECTED :: edge_opt = 0._xp ! meant to redistribute the points in z
+ REAL(xp), PUBLIC, PROTECTED :: major_R = 1._xp ! major radius
+ REAL(xp), PUBLIC, PROTECTED :: major_Z = 0._xp ! vertical elevation
+ REAL(xp), PUBLIC, PROTECTED :: xpdx_pm_geom = 0._xp ! amplitude mag. eq. pressure grad.
+ REAL(xp), PUBLIC, PROTECTED :: C_y = 0._xp ! defines y coordinate : Cy (q theta - phi)
+ REAL(xp), PUBLIC, PROTECTED :: C_xy = 1._xp ! defines x coordinate : B = Cxy Vx x Vy
! Geometrical auxiliary variables
LOGICAL, PUBLIC, PROTECTED :: SHEARED = .false. ! flag for shear magn. geom or not
! Curvature
- REAL(dp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z, odd/even p
+ REAL(xp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z, odd/even p
! Jacobian
- REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Jacobian ! dimensions: z, odd/even p
- COMPLEX(dp), PUBLIC, PROTECTED :: iInt_Jacobian ! Inverse integrated Jacobian
+ REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Jacobian ! dimensions: z, odd/even p
+ COMPLEX(xp), PUBLIC, PROTECTED :: iInt_Jacobian ! Inverse integrated Jacobian
! Metric
- REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gxx, gxy, gxz, gyy, gyz, gzz ! dimensions: z, odd/even p
- REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dxdr, dxdZ, Rc, phic, Zc
+ REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gxx, gxy, gxz, gyy, gyz, gzz ! dimensions: z, odd/even p
+ REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dxdr, dxdZ, Rc, phic, Zc
! derivatives of magnetic field strength
- REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dBdx, dBdy, dBdz, dlnBdz
+ REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dBdx, dBdy, dBdz, dlnBdz
! Relative magnetic field strength
- REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB
+ REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB
! Relative strength of major radius
- REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ
+ REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ
! Some geometrical coefficients
- REAL(dp), PUBLIC, DIMENSION(:,:) , ALLOCATABLE :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
+ REAL(xp), PUBLIC, DIMENSION(:,:) , ALLOCATABLE :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
! Array to map the index of mode (kx,ky,-pi) to (kx+2pi*s*ky,ky,pi) for sheared periodic boundary condition
INTEGER, PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ikx_zBC_L, ikx_zBC_R
! Geometric factor in front of the parallel phi derivative (not implemented)
- ! REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Gamma_phipar
+ ! REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Gamma_phipar
! pb_phase, for parallel boundary phase, contains the factor that occurs when taking into account
! that q0 is defined in the middle of the fluxtube whereas the radial position spans in [0,Lx)
! This shift introduces a (-1)^(Nexc*iky) phase change that is included in GENE
- COMPLEX(dp), PUBLIC, DIMENSION(:), ALLOCATABLE :: pb_phase_L, pb_phase_R
+ COMPLEX(xp), PUBLIC, DIMENSION(:), ALLOCATABLE :: pb_phase_L, pb_phase_R
! Functions
PUBLIC :: geometry_readinputs, geometry_outputinputs,&
@@ -78,7 +78,7 @@ CONTAINS
kappa, s_kappa,delta, s_delta, zeta, s_zeta,& ! For miller
parallel_bc, shift_y, Npol
READ(lu_in,geometry)
- IF(shear .NE. 0._dp) SHEARED = .true.
+ IF(shear .NE. 0._xp) SHEARED = .true.
SELECT CASE(parallel_bc)
CASE ('dirichlet')
CASE ('periodic')
@@ -100,9 +100,9 @@ CONTAINS
USE calculus, ONLY: simpson_rule_z
! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient
implicit none
- REAL(dp) :: kx,ky
- COMPLEX(dp), DIMENSION(local_nz) :: integrant
- real(dp) :: G1,G2,G3,Cx,Cy
+ REAL(xp) :: kx,ky
+ COMPLEX(xp), DIMENSION(local_nz) :: integrant
+ real(xp) :: G1,G2,G3,Cx,Cy
INTEGER :: eo,iz,iky,ikx
! Allocate arrays
@@ -120,12 +120,12 @@ CONTAINS
CALL speak('Z-pinch geometry')
call eval_zpinch_geometry
SHEARED = .FALSE.
- shear = 0._dp
+ shear = 0._xp
CASE('miller')
CALL speak('Miller geometry')
call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta)
call get_miller(eps,major_R,major_Z,q0,shear,Npol,alpha_MHD,edge_opt,&
- C_y,C_xy,dpdx_pm_geom,gxx,gyy,gzz,gxy,gxz,gyz,&
+ C_y,C_xy,xpdx_pm_geom,gxx,gyy,gzz,gxy,gxz,gyz,&
dBdx,dBdy,hatB,jacobian,dBdz,hatR,hatZ,dxdR,dxdZ,&
Ckxky,gradz_coeff)
CASE DEFAULT
@@ -154,7 +154,7 @@ CONTAINS
ENDDO
ENDDO
! coefficient in the front of parallel derivative
- gradz_coeff(iz,eo) = 1._dp /(jacobian(iz,eo)*hatB(iz,eo))
+ gradz_coeff(iz,eo) = 1._xp /(jacobian(iz,eo)*hatB(iz,eo))
! d/dz(ln B) to correspond to the formulation in paper 2023
dlnBdz(iz,eo) = dBdz(iz,eo)/hatB(iz,eo)
! Geometric factor in front to the maxwellian dzphi term (not implemented)
@@ -168,7 +168,7 @@ CONTAINS
! Compute the inverse z integrated Jacobian (useful for flux averaging)
integrant = Jacobian((1+ngz/2):(local_nz+ngz/2),ieven) ! Convert into complex array
CALL simpson_rule_z(local_nz,zweights_SR,integrant,iInt_Jacobian)
- iInt_Jacobian = 1._dp/iInt_Jacobian ! reverse it
+ iInt_Jacobian = 1._xp/iInt_Jacobian ! reverse it
END SUBROUTINE eval_magnetic_geometry
!
!--------------------------------------------------------------------------------
@@ -178,27 +178,27 @@ CONTAINS
USE grid, ONLY : local_nz,Ngz,zarray,nzgrid
! evaluate s-alpha geometry model
implicit none
- REAL(dp) :: z
+ REAL(xp) :: z
INTEGER :: iz, eo
- alpha_MHD = 0._dp
+ alpha_MHD = 0._xp
DO eo = 1,nzgrid
DO iz = 1,local_nz+Ngz
z = zarray(iz,eo)
! metric
- gxx(iz,eo) = 1._dp
+ gxx(iz,eo) = 1._xp
gxy(iz,eo) = shear*z - alpha_MHD*SIN(z)
- gxz(iz,eo) = 0._dp
- gyy(iz,eo) = 1._dp + (shear*z - alpha_MHD*SIN(z))**2
- gyz(iz,eo) = 1._dp/eps
- gzz(iz,eo) = 1._dp/eps**2
+ gxz(iz,eo) = 0._xp
+ gyy(iz,eo) = 1._xp + (shear*z - alpha_MHD*SIN(z))**2
+ gyz(iz,eo) = 1._xp/eps
+ gzz(iz,eo) = 1._xp/eps**2
dxdR(iz,eo)= COS(z)
dxdZ(iz,eo)= SIN(z)
! Poloidal plane coordinates
- hatR(iz,eo) = 1._dp + eps*COS(z)
- hatZ(iz,eo) = 1._dp + eps*SIN(z)
+ hatR(iz,eo) = 1._xp + eps*COS(z)
+ hatZ(iz,eo) = 1._xp + eps*SIN(z)
! toroidal coordinates
Rc (iz,eo) = hatR(iz,eo)
@@ -206,18 +206,18 @@ CONTAINS
Zc (iz,eo) = hatZ(iz,eo)
! Relative strengh of modulus of B
- hatB(iz,eo) = 1._dp/(1._dp + eps*COS(z))
+ hatB(iz,eo) = 1._xp/(1._xp + eps*COS(z))
! Jacobian
Jacobian(iz,eo) = q0/hatB(iz,eo)
! Derivative of the magnetic field strenght
dBdx(iz,eo) = -COS(z)*hatB(iz,eo)**2 ! LB = 1
- dBdy(iz,eo) = 0._dp
+ dBdy(iz,eo) = 0._xp
dBdz(iz,eo) = eps*SIN(z)*hatB(iz,eo)**2
! Curvature factor
- C_xy = 1._dp
+ C_xy = 1._xp
ENDDO
ENDDO
@@ -229,27 +229,27 @@ CONTAINS
SUBROUTINE eval_zpinch_geometry
USE grid, ONLY : local_nz,Ngz,zarray,nzgrid
implicit none
- REAL(dp) :: z
+ REAL(xp) :: z
INTEGER :: iz, eo
- alpha_MHD = 0._dp
+ alpha_MHD = 0._xp
DO eo = 1,nzgrid
DO iz = 1,local_nz+Ngz
z = zarray(iz,eo)
! metric
- gxx(iz,eo) = 1._dp
- gxy(iz,eo) = 0._dp
- gxz(iz,eo) = 0._dp
- gyy(iz,eo) = 1._dp ! 1/R but R is the normalization length
- gyz(iz,eo) = 0._dp
- gzz(iz,eo) = 1._dp
+ gxx(iz,eo) = 1._xp
+ gxy(iz,eo) = 0._xp
+ gxz(iz,eo) = 0._xp
+ gyy(iz,eo) = 1._xp ! 1/R but R is the normalization length
+ gyz(iz,eo) = 0._xp
+ gzz(iz,eo) = 1._xp
dxdR(iz,eo)= COS(z)
dxdZ(iz,eo)= SIN(z)
! Relative strengh of radius
- hatR(iz,eo) = 1._dp ! R but R is the normalization length
- hatZ(iz,eo) = 1._dp
+ hatR(iz,eo) = 1._xp ! R but R is the normalization length
+ hatZ(iz,eo) = 1._xp
! toroidal coordinates
Rc (iz,eo) = hatR(iz,eo)
@@ -257,21 +257,21 @@ CONTAINS
Zc (iz,eo) = hatZ(iz,eo)
! Jacobian
- Jacobian(iz,eo) = 1._dp ! R but R is the normalization length
+ Jacobian(iz,eo) = 1._xp ! R but R is the normalization length
! Relative strengh of modulus of B
- hatB (iz,eo) = 1._dp
+ hatB (iz,eo) = 1._xp
! Derivative of the magnetic field strenght
dBdx(iz,eo) = -hatB(iz,eo) ! LB = 1
- dBdy(iz,eo) = 0._dp
- dBdz(iz,eo) = 0._dp ! Gene put a factor hatB or 1/hatR in this
+ dBdy(iz,eo) = 0._xp
+ dBdz(iz,eo) = 0._xp ! Gene put a factor hatB or 1/hatR in this
ENDDO
ENDDO
! Curvature factor
- C_xy = 1._dp
+ C_xy = 1._xp
END SUBROUTINE eval_zpinch_geometry
!
!--------------------------------------------------------------------------------
@@ -280,25 +280,25 @@ CONTAINS
USE grid, ONLY : local_nz,Ngz,zarray, nzgrid
! evaluate 1D perp geometry model
implicit none
- REAL(dp) :: z
+ REAL(xp) :: z
INTEGER :: iz, eo
DO eo = 1,nzgrid
DO iz = 1,local_nz+Ngz
z = zarray(iz,eo)
! metric
- gxx(iz,eo) = 1._dp
- gxy(iz,eo) = 0._dp
- gyy(iz,eo) = 1._dp
+ gxx(iz,eo) = 1._xp
+ gxy(iz,eo) = 0._xp
+ gyy(iz,eo) = 1._xp
! Relative strengh of radius
- hatR(iz,eo) = 1._dp
+ hatR(iz,eo) = 1._xp
! Jacobian
- Jacobian(iz,eo) = 1._dp
+ Jacobian(iz,eo) = 1._xp
! Relative strengh of modulus of B
- hatB(iz,eo) = 1._dp
+ hatB(iz,eo) = 1._xp
ENDDO
ENDDO
@@ -313,7 +313,7 @@ CONTAINS
local_nky_offset
USE prec_const, ONLY: imagu, pi
IMPLICIT NONE
- ! REAL(dp) :: shift
+ ! REAL(xp) :: shift
INTEGER :: ikx,iky, mn_y
ALLOCATE(ikx_zBC_L(local_nky,total_nkx))
ALLOCATE(ikx_zBC_R(local_nky,total_nkx))
@@ -330,8 +330,8 @@ CONTAINS
ikx_zBC_L(iky,ikx) = ikx ! connect to itself per default
ikx_zBC_R(iky,ikx) = ikx
ENDDO
- pb_phase_L(iky) = 1._dp ! no phase change per default
- pb_phase_R(iky) = 1._dp
+ pb_phase_L(iky) = 1._xp ! no phase change per default
+ pb_phase_R(iky) = 1._xp
ENDDO
! Parallel boundary are not trivial for sheared case and if
! the user does not ask explicitly for shearless bc
@@ -343,7 +343,7 @@ CONTAINS
! get the real mode number (iky starts at 1 and is shifted from paral)
mn_y = iky-1+local_nky_offset
! Formula for the shift due to shear after Npol turns
- ! shift = 2._dp*PI*shear*kyarray(iky)*Npol
+ ! shift = 2._xp*PI*shear*kyarray(iky)*Npol
DO ikx = 1,total_nkx
! Usual formula for shifting indices using that dkx = 2pi*shear*dky/Nexc
ikx_zBC_L(iky,ikx) = ikx-mn_y*Nexc
@@ -362,7 +362,7 @@ CONTAINS
ENDDO
! phase present in GENE from a shift of the x origin by Lx/2 (useless?)
! We also put the user defined shift in the y direction (see Volcokas et al. 2022)
- pb_phase_L(iky) = (-1._dp)**(Nexc*mn_y)*EXP(imagu*REAL(mn_y,dp)*2._dp*pi*shift_y)
+ pb_phase_L(iky) = (-1._xp)**(Nexc*mn_y)*EXP(imagu*REAL(mn_y,xp)*2._xp*pi*shift_y)
ENDDO
ENDIF
! Option for disconnecting every modes, viz. connecting all boundary to 0
@@ -373,7 +373,7 @@ CONTAINS
! get the real mode number (iky starts at 1 and is shifted from paral)
mn_y = iky-1+local_nky_offset
! Formula for the shift due to shear after Npol
- ! shift = 2._dp*PI*shear*kyarray(iky)*Npol
+ ! shift = 2._xp*PI*shear*kyarray(iky)*Npol
DO ikx = 1,total_nkx
! Usual formula for shifting indices
ikx_zBC_R(iky,ikx) = ikx+mn_y*Nexc
@@ -393,7 +393,7 @@ CONTAINS
ENDDO
! phase present in GENE from a shift ofthe x origin by Lx/2 (useless?)
! We also put the user defined shift in the y direction (see Volcokas et al. 2022)
- pb_phase_R(iky) = (-1._dp)**(Nexc*mn_y)*EXP(-imagu*REAL(mn_y,dp)*2._dp*pi*shift_y)
+ pb_phase_R(iky) = (-1._xp)**(Nexc*mn_y)*EXP(-imagu*REAL(mn_y,xp)*2._xp*pi*shift_y)
ENDDO
ENDIF
! Option for disconnecting every modes, viz. connecting all boundary to 0
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index c578ad1b..9a2279b8 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -1,6 +1,6 @@
module ghosts
USE mpi
-USE prec_const, ONLY: dp
+USE prec_const, ONLY: xp
IMPLICIT NONE
INTEGER :: status(MPI_STATUS_SIZE), source, dest, count, ipg
@@ -28,7 +28,7 @@ SUBROUTINE update_ghosts_EM
IMPLICIT NONE
IF(total_nz .GT. 1) THEN
CALL update_ghosts_z_3D(phi)
- IF(beta .GT. 0._dp) &
+ IF(beta .GT. 0._xp) &
CALL update_ghosts_z_3D(psi)
ENDIF
END SUBROUTINE update_ghosts_EM
@@ -102,7 +102,7 @@ SUBROUTINE update_ghosts_z_mom
ngp,ngj,ngz
IMPLICIT NONE
!! buffer for data exchanges
- COMPLEX(dp),DIMENSION(local_na,local_np+ngp,local_nj+ngj,local_nky,local_nkx,-Ngz/2:Ngz/2) :: buff_pjxy_zBC
+ COMPLEX(xp),DIMENSION(local_na,local_np+ngp,local_nj+ngj,local_nky,local_nkx,-Ngz/2:Ngz/2) :: buff_pjxy_zBC
INTEGER :: ikxBC_L, ikxBC_R, ikx, iky, first, last, ig, ierr
first = 1 + ngz/2
last = local_nz + ngz/2
@@ -139,7 +139,7 @@ SUBROUTINE update_ghosts_z_mom
ENDDO
ELSE
DO ig=1,ngz/2
- moments(:,:,:,iky,ikx,first-ig,updatetlevel) = 0._dp
+ moments(:,:,:,iky,ikx,first-ig,updatetlevel) = 0._xp
ENDDO
ENDIF
ikxBC_R = ikx_zBC_R(iky,ikx);
@@ -151,7 +151,7 @@ SUBROUTINE update_ghosts_z_mom
ENDDO
ELSE
DO ig=1,ngz/2
- moments(:,:,:,iky,ikx,last+ig,updatetlevel) = 0._dp
+ moments(:,:,:,iky,ikx,last+ig,updatetlevel) = 0._xp
ENDDO
ENDIF
ENDDO
@@ -164,8 +164,8 @@ SUBROUTINE update_ghosts_z_3D(field)
USE grid, ONLY: local_nky,local_nkx,local_nz,ngz
IMPLICIT NONE
!! buffer for data exchanges
- COMPLEX(dp),DIMENSION(local_nky,local_nkx,-ngz/2:ngz/2) :: buff_xy_zBC
- COMPLEX(dp), INTENT(INOUT) :: field(local_nky,local_nkx,local_nz+ngz)
+ COMPLEX(xp),DIMENSION(local_nky,local_nkx,-ngz/2:ngz/2) :: buff_xy_zBC
+ COMPLEX(xp), INTENT(INOUT) :: field(local_nky,local_nkx,local_nz+ngz)
INTEGER :: ikxBC_L, ikxBC_R, ikx, iky, first, last, ig, ierr
first = 1 + ngz/2
last = local_nz + ngz/2
@@ -200,7 +200,7 @@ SUBROUTINE update_ghosts_z_3D(field)
ENDDO
ELSE
DO ig = 1,ngz/2
- field(iky,ikx,first-ig) = 0._dp
+ field(iky,ikx,first-ig) = 0._xp
ENDDO
ENDIF
ikxBC_R = ikx_zBC_R(iky,ikx);
@@ -211,7 +211,7 @@ SUBROUTINE update_ghosts_z_3D(field)
ENDDO
ELSE
DO ig = 1,ngz/2
- field(iky,ikx,last+ig) = 0._dp
+ field(iky,ikx,last+ig) = 0._xp
ENDDO
ENDIF
ENDDO
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 95b9c588..671f015f 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -13,23 +13,23 @@ MODULE grid
INTEGER, PUBLIC, PROTECTED :: maxj = 1 ! The maximal Laguerre-moment
INTEGER, PUBLIC, PROTECTED :: dmax = 1 ! The maximal full GF set of i-moments v^dmax
INTEGER, PUBLIC, PROTECTED :: Nx = 4 ! Number of total internal grid points in x
- REAL(dp), PUBLIC, PROTECTED :: Lx = 120_dp ! horizontal length of the spatial box
+ REAL(xp), PUBLIC, PROTECTED :: Lx = 120_xp ! horizontal length of the spatial box
INTEGER, PUBLIC, PROTECTED :: Nexc = 1 ! factor to increase Lx when shear>0 (Lx = Nexc/kymin/shear)
INTEGER, PUBLIC, PROTECTED :: Ny = 4 ! Number of total internal grid points in y
- REAL(dp), PUBLIC, PROTECTED :: Ly = 120_dp ! vertical length of the spatial box
+ REAL(xp), PUBLIC, PROTECTED :: Ly = 120_xp ! vertical length of the spatial box
INTEGER, PUBLIC, PROTECTED :: Nz = 4 ! Number of total perpendicular planes
INTEGER, PUBLIC, PROTECTED :: Odz = 4 ! order of z interp and derivative schemes
INTEGER, PUBLIC, PROTECTED :: Nkx ! Number of total internal grid points in kx
INTEGER, PUBLIC, PROTECTED :: Nky ! Number of total internal grid points in ky
- REAL(dp), PUBLIC, PROTECTED :: kpar = 0_dp ! parallel wave vector component
+ REAL(xp), PUBLIC, PROTECTED :: kpar = 0_xp ! parallel wave vector component
! Grid arrays
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: parray, parray_full
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: jarray, jarray_full
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kxarray, kxarray_full
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kyarray, kyarray_full
- REAL(dp), DIMENSION(:,:), ALLOCATABLE, PUBLIC,PROTECTED :: zarray
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: zarray_full
- REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE, PUBLIC,PROTECTED :: kparray !kperp
+ REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kxarray, kxarray_full
+ REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kyarray, kyarray_full
+ REAL(xp), DIMENSION(:,:), ALLOCATABLE, PUBLIC,PROTECTED :: zarray
+ REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: zarray_full
+ REAL(xp), DIMENSION(:,:,:,:), ALLOCATABLE, PUBLIC,PROTECTED :: kparray !kperp
! Indexation variables
INTEGER, PUBLIC, PROTECTED :: ias ,iae ! species index
INTEGER, PUBLIC, PROTECTED :: ips ,ipe ! Hermite
@@ -68,18 +68,18 @@ MODULE grid
integer(C_INTPTR_T), PUBLIC,PROTECTED :: local_nkx_ptr, local_nky_ptr
integer(C_INTPTR_T), PUBLIC,PROTECTED :: local_nkx_ptr_offset, local_nky_ptr_offset
! Grid spacing and limits
- REAL(dp), PUBLIC, PROTECTED :: deltap, deltaz, inv_deltaz
- REAL(dp), PUBLIC, PROTECTED :: deltakx, deltaky, kx_max, ky_max, kx_min, ky_min!, kp_max
+ REAL(xp), PUBLIC, PROTECTED :: deltap, deltaz, inv_deltaz
+ REAL(xp), PUBLIC, PROTECTED :: deltakx, deltaky, kx_max, ky_max, kx_min, ky_min!, kp_max
INTEGER , PUBLIC, PROTECTED :: local_pmin, local_pmax
INTEGER , PUBLIC, PROTECTED :: local_jmin, local_jmax
- REAL(dp), PUBLIC, PROTECTED :: local_kymin, local_kymax
- REAL(dp), PUBLIC, PROTECTED :: local_kxmin, local_kxmax
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: local_zmin, local_zmax
+ REAL(xp), PUBLIC, PROTECTED :: local_kymin, local_kymax
+ REAL(xp), PUBLIC, PROTECTED :: local_kxmin, local_kxmax
+ REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: local_zmin, local_zmax
! local z weights for computing simpson rule
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: zweights_SR
+ REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: zweights_SR
! Numerical diffusion scaling
- REAL(dp), PUBLIC, PROTECTED :: diff_p_coeff, diff_j_coeff
- REAL(dp), PUBLIC, PROTECTED :: diff_kx_coeff, diff_ky_coeff, diff_dz_coeff
+ REAL(xp), PUBLIC, PROTECTED :: diff_p_coeff, diff_j_coeff
+ REAL(xp), PUBLIC, PROTECTED :: diff_kx_coeff, diff_ky_coeff, diff_dz_coeff
LOGICAL, PUBLIC, PROTECTED :: SG = .true.! shifted grid flag
! Array to know the distribution of data among all processes (for MPI comm)
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: counts_total_nkx, counts_nky, counts_nz
@@ -95,14 +95,14 @@ MODULE grid
LOGICAL, PUBLIC, PROTECTED :: CONTAINSp0, CONTAINSp1, CONTAINSp2, CONTAINSp3
LOGICAL, PUBLIC, PROTECTED :: SOLVE_POISSON, SOLVE_AMPERE
! Usefull inverse numbers
- REAL(dp), PUBLIC, PROTECTED :: inv_Nx, inv_Ny, inv_Nz
+ REAL(xp), PUBLIC, PROTECTED :: inv_Nx, inv_Ny, inv_Nz
! For Orszag filter
- REAL(dp), PUBLIC, PROTECTED :: two_third_kxmax
- REAL(dp), PUBLIC, PROTECTED :: two_third_kymax
- REAL(dp), PUBLIC, PROTECTED :: two_third_kpmax
+ REAL(xp), PUBLIC, PROTECTED :: two_third_kxmax
+ REAL(xp), PUBLIC, PROTECTED :: two_third_kymax
+ REAL(xp), PUBLIC, PROTECTED :: two_third_kpmax
! 1D Antialiasing arrays (2/3 rule)
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: AA_x
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: AA_y
+ REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: AA_x
+ REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: AA_y
! Public Functions
PUBLIC :: init_1Dgrid_distr
@@ -111,7 +111,7 @@ MODULE grid
PUBLIC :: bar
! Precomputations
- real(dp), PUBLIC, PROTECTED :: pmax_dp, jmax_dp
+ real(xp), PUBLIC, PROTECTED :: pmax_xp, jmax_xp
CONTAINS
@@ -129,8 +129,8 @@ CONTAINS
! i.e. : all moments N_a^pj s.t. p+2j<=d are simulated (see GF closure)
dmax = min(pmax,2*jmax+1)
! Usefull precomputations
- inv_Nx = 1._dp/REAL(Nx,dp)
- inv_Ny = 1._dp/REAL(Ny,dp)
+ inv_Nx = 1._xp/REAL(Nx,xp)
+ inv_Ny = 1._xp/REAL(Ny,xp)
END SUBROUTINE grid_readinputs
!!!! GRID REPRESENTATION
! We define the grids that contain ghosts (p,j,z) with indexing from 1 to Nlocal + nghost
@@ -143,7 +143,7 @@ CONTAINS
! 1 2 3 4
SUBROUTINE set_grids(shear,Npol,LINEARITY,N_HD,EM,Na)
USE fourier, ONLY: init_grid_distr_and_plans
- REAL(dp), INTENT(IN) :: shear
+ REAL(xp), INTENT(IN) :: shear
INTEGER, INTENT(IN) :: Npol
CHARACTER(len=*), INTENT(IN) :: LINEARITY
INTEGER, INTENT(IN) :: N_HD
@@ -252,8 +252,8 @@ CONTAINS
IF(CONTAINSp0 .AND. .NOT. (CONTAINSp2))&
WRITE(*,*) "Warning : distribution along p may not work with DGGK"
! Precomputations
- pmax_dp = real(pmax,dp)
- diff_p_coeff = pmax_dp*(1._dp/pmax_dp)**6
+ pmax_xp = real(pmax,xp)
+ diff_p_coeff = pmax_xp*(1._xp/pmax_xp)**6
! Overwrite SOLVE_AMPERE flag if beta is zero
IF(.NOT. EM) THEN
SOLVE_AMPERE = .FALSE.
@@ -284,8 +284,8 @@ CONTAINS
local_jmax = jarray(local_nj+ngj/2)
local_jmin = jarray(1+ngj/2)
! Precomputations
- jmax_dp = real(jmax,dp)
- diff_j_coeff = jmax_dp*(1._dp/jmax_dp)**6
+ jmax_xp = real(jmax,xp)
+ diff_j_coeff = jmax_xp*(1._xp/jmax_xp)**6
! j=0 and j=1 indices
DO ij = 1,local_nj+ngj
IF(jarray(ij) .EQ. 0) ij0 = ij
@@ -303,11 +303,11 @@ CONTAINS
total_nky = Nky
! Grid spacings
IF (Ny .EQ. 1) THEN
- deltaky = 2._dp*PI/Ly
+ deltaky = 2._xp*PI/Ly
ky_max = deltaky
ky_min = deltaky
ELSE
- deltaky = 2._dp*PI/Ly
+ deltaky = 2._xp*PI/Ly
ky_max = (Nky-1)*deltaky
ky_min = deltaky
ENDIF
@@ -315,14 +315,14 @@ CONTAINS
! Build the full grids on process 0 to diagnose it without comm
ALLOCATE(kyarray_full(Nky))
DO iky = 1,Nky
- kyarray_full(iky) = REAL(iky-1,dp) * deltaky
+ kyarray_full(iky) = REAL(iky-1,xp) * deltaky
END DO
ikys = local_nky_ptr_offset + 1
ikye = ikys + local_nky_ptr - 1
local_nky = ikye - ikys + 1
local_nky_offset = local_nky_ptr_offset
ALLOCATE(kyarray(local_nky))
- local_kymax = 0._dp
+ local_kymax = 0._xp
! Creating a grid ordered as dk*(0 1 2 3)
! We loop over the natural iky numbers (|1 2 3||4 5 6||... Nky|)
DO iky = 1,local_nky
@@ -351,36 +351,36 @@ CONTAINS
ENDIF
END DO
! Orszag 2/3 filter
- two_third_kymax = 2._dp/3._dp*deltaky*(Nky-1)
+ two_third_kymax = 2._xp/3._xp*deltaky*(Nky-1)
ALLOCATE(AA_y(local_nky))
DO iky = 1,local_nky
IF ( (kyarray(iky) .LT. two_third_kymax) .OR. (LINEARITY .EQ. 'linear')) THEN
- AA_y(iky) = 1._dp;
+ AA_y(iky) = 1._xp;
ELSE
- AA_y(iky) = 0._dp;
+ AA_y(iky) = 0._xp;
ENDIF
END DO
! For hyperdiffusion
IF(LINEARITY.EQ.'linear') THEN
- diff_ky_coeff= (1._dp/ky_max)**N_HD
+ diff_ky_coeff= (1._xp/ky_max)**N_HD
ELSE
- diff_ky_coeff= (1._dp/two_third_kymax)**N_HD
+ diff_ky_coeff= (1._xp/two_third_kymax)**N_HD
ENDIF
END SUBROUTINE set_kygrid
SUBROUTINE set_kxgrid(shear,Npol,LINEARITY,N_HD)
USE prec_const
IMPLICIT NONE
- REAL(dp), INTENT(IN) :: shear
+ REAL(xp), INTENT(IN) :: shear
INTEGER, INTENT(IN) :: Npol
CHARACTER(len=*), INTENT(IN) ::LINEARITY
INTEGER, INTENT(IN) :: N_HD
INTEGER :: ikx
- REAL(dp):: Lx_adapted
+ REAL(xp):: Lx_adapted
IF(shear .GT. 0) THEN
IF(my_id.EQ.0) write(*,*) 'Magnetic shear detected: set up sheared kx grid..'
! mininal size of box in x to respect dkx = 2pi shear dky
- Lx_adapted = Ly/(2._dp*pi*shear*Npol)
+ Lx_adapted = Ly/(2._xp*pi*shear*Npol)
! Put Nexc to 0 so that it is computed from a target value Lx
IF(Nexc .EQ. 0) THEN
Nexc = CEILING(0.9 * Lx/Lx_adapted)
@@ -401,27 +401,27 @@ CONTAINS
ALLOCATE(kxarray(local_nkx))
ALLOCATE(kxarray_full(total_nkx))
IF (Nx .EQ. 1) THEN
- deltakx = 1._dp
- kxarray(1) = 0._dp
+ deltakx = 1._xp
+ kxarray(1) = 0._xp
ikx0 = 1
contains_kx0 = .true.
- kx_max = 0._dp
+ kx_max = 0._xp
ikx_max = 1
- kx_min = 0._dp
- kxarray_full(1) = 0._dp
- local_kxmax = 0._dp
+ kx_min = 0._xp
+ kxarray_full(1) = 0._xp
+ local_kxmax = 0._xp
ELSE ! Build apprpopriate grid
- deltakx = 2._dp*PI/Lx
+ deltakx = 2._xp*PI/Lx
IF(MODULO(total_nkx,2) .EQ. 0) THEN ! Even number of kx (-2 -1 0 1 2 3)
kx_max = (total_nkx/2)*deltakx
kx_min = -kx_max+deltakx
! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
DO ikx = 1,total_nkx
- kxarray_full(ikx) = deltakx*REAL(MODULO(ikx-1,total_nkx/2)-(total_nkx/2)*FLOOR(2.*real(ikx-1)/real(total_nkx)),dp)
+ kxarray_full(ikx) = deltakx*REAL(MODULO(ikx-1,total_nkx/2)-(total_nkx/2)*FLOOR(2.*real(ikx-1)/real(total_nkx)),xp)
IF (ikx .EQ. total_nkx/2+1) kxarray_full(ikx) = -kxarray_full(ikx)
END DO
! Set local grid (not parallelized so same as full one)
- local_kxmax = 0._dp
+ local_kxmax = 0._xp
DO ikx = 1,local_nkx
kxarray(ikx) = kxarray_full(ikx+local_nkx_offset)
! Finding kx=0
@@ -440,22 +440,22 @@ CONTAINS
ENDIF
ENDIF
! Orszag 2/3 filter
- two_third_kxmax = 2._dp/3._dp*kx_max;
+ two_third_kxmax = 2._xp/3._xp*kx_max;
! Antialiasing filter
ALLOCATE(AA_x(local_nkx))
DO ikx = 1,local_nkx
IF ( ((kxarray(ikx) .GT. -two_third_kxmax) .AND. &
(kxarray(ikx) .LT. two_third_kxmax)) .OR. (LINEARITY .EQ. 'linear')) THEN
- AA_x(ikx) = 1._dp;
+ AA_x(ikx) = 1._xp;
ELSE
- AA_x(ikx) = 0._dp;
+ AA_x(ikx) = 0._xp;
ENDIF
END DO
! For hyperdiffusion
IF(LINEARITY.EQ.'linear') THEN
- diff_kx_coeff= (1._dp/kx_max)**N_HD
+ diff_kx_coeff= (1._xp/kx_max)**N_HD
ELSE
- diff_kx_coeff= (1._dp/two_third_kxmax)**N_HD
+ diff_kx_coeff= (1._xp/two_third_kxmax)**N_HD
ENDIF
END SUBROUTINE set_kxgrid
@@ -463,25 +463,25 @@ CONTAINS
USE prec_const
USE parallel, ONLY: num_procs_z, rank_z
IMPLICIT NONE
- REAL(dp):: grid_shift, Lz, zmax, zmin
+ REAL(xp):: grid_shift, Lz, zmax, zmin
INTEGER :: istart, iend, in, Npol, iz, ig, eo, iglob
total_nz = Nz
! Length of the flux tube (in ballooning angle)
- Lz = 2._dp*pi*REAL(Npol,dp)
+ Lz = 2._xp*pi*REAL(Npol,xp)
! Z stepping (#interval = #points since periodic)
- deltaz = Lz/REAL(Nz,dp)
- inv_deltaz = 1._dp/deltaz
+ deltaz = Lz/REAL(Nz,xp)
+ inv_deltaz = 1._xp/deltaz
! Parallel hyperdiffusion coefficient
- diff_dz_coeff = (deltaz/2._dp)**4 ! adaptive fourth derivative (~GENE)
+ diff_dz_coeff = (deltaz/2._xp)**4 ! adaptive fourth derivative (~GENE)
IF (SG) THEN
CALL speak('--2 staggered z grids--')
- grid_shift = deltaz/2._dp
+ grid_shift = deltaz/2._xp
! indices for even p and odd p grids (used in kernel, jacobian, gij etc.)
ieven = 1
iodd = 2
nzgrid = 2
ELSE
- grid_shift = 0._dp
+ grid_shift = 0._xp
ieven = 1
iodd = 1
nzgrid = 1
@@ -491,7 +491,7 @@ CONTAINS
IF (Nz .EQ. 1) Npol = 0
zmax = 0; zmin = 0;
DO iz = 1,total_nz ! z in [-pi pi-dz] x Npol
- zarray_full(iz) = REAL(iz-1,dp)*deltaz - Lz/2._dp
+ zarray_full(iz) = REAL(iz-1,xp)*deltaz - Lz/2._xp
IF(zarray_full(iz) .GT. zmax) zmax = zarray_full(iz)
IF(zarray_full(iz) .LT. zmin) zmin = zarray_full(iz)
END DO
@@ -527,7 +527,7 @@ CONTAINS
!! interior point loop
DO iz = 1,local_nz
DO eo = 1,nzgrid
- zarray(iz+ngz/2,eo) = zarray_full(iz+local_nz_offset) + REAL(eo-1,dp)*grid_shift
+ zarray(iz+ngz/2,eo) = zarray_full(iz+local_nz_offset) + REAL(eo-1,xp)*grid_shift
ENDDO
ENDDO
CALL allocate_array(local_zmax,1,nzgrid)
@@ -539,8 +539,8 @@ CONTAINS
! Fill the ghosts
! Continue angles
! DO ig = 1,ngz/2
- ! zarray(ig,eo) = local_zmin(eo)-REAL(ngz/2-(ig-1),dp)*deltaz
- ! zarray(local_nz+ngz/2+ig,eo) = local_zmax(eo)+REAL(ig,dp)*deltaz
+ ! zarray(ig,eo) = local_zmin(eo)-REAL(ngz/2-(ig-1),xp)*deltaz
+ ! zarray(local_nz+ngz/2+ig,eo) = local_zmax(eo)+REAL(ig,xp)*deltaz
! ENDDO
! Periodic z \in (-pi pi-dz)
DO ig = 1,ngz/2 ! first ghost cells
@@ -557,30 +557,30 @@ CONTAINS
ENDDO
! Set up the flags to know if the process contains the tip and/or the tail
! of the z domain (important for z-boundary condition)
- IF(abs(local_zmin(eo) - (zmin+REAL(eo-1,dp)*grid_shift)) .LT. EPSILON(zmin)) &
+ IF(abs(local_zmin(eo) - (zmin+REAL(eo-1,xp)*grid_shift)) .LT. EPSILON(zmin)) &
contains_zmin = .TRUE.
- IF(abs(local_zmax(eo) - (zmax+REAL(eo-1,dp)*grid_shift)) .LT. EPSILON(zmax)) &
+ IF(abs(local_zmax(eo) - (zmax+REAL(eo-1,xp)*grid_shift)) .LT. EPSILON(zmax)) &
contains_zmax = .TRUE.
ENDDO
! local weights for Simpson rule
ALLOCATE(zweights_SR(local_nz))
IF(total_nz .EQ. 1) THEN
- zweights_SR = 1._dp
+ zweights_SR = 1._xp
ELSE
DO iz = 1,local_nz
IF(MODULO(iz+local_nz_offset,2) .EQ. 1) THEN ! odd iz
- zweights_SR(iz) = onethird*deltaz*2._dp
+ zweights_SR(iz) = onethird*deltaz*2._xp
ELSE ! even iz
- zweights_SR(iz) = onethird*deltaz*4._dp
+ zweights_SR(iz) = onethird*deltaz*4._xp
ENDIF
ENDDO
ENDIF
END SUBROUTINE set_zgrid
SUBROUTINE set_kparray(gxx, gxy, gyy,hatB)
- REAL(dp), DIMENSION(local_nz+ngz,nzgrid), INTENT(IN) :: gxx,gxy,gyy,hatB
+ REAL(xp), DIMENSION(local_nz+ngz,nzgrid), INTENT(IN) :: gxx,gxy,gyy,hatB
INTEGER :: eo,iz,iky,ikx
- REAL(dp) :: kx, ky
+ REAL(xp) :: kx, ky
CALL allocate_array( kparray, 1,local_nky, 1,local_nkx, 1,local_nz+ngz, 1,nzgrid)
DO eo = 1,nzgrid
DO iz = 1,local_nz+ngz
@@ -591,12 +591,12 @@ CONTAINS
! there is a factor 1/B from the normalization; important to match GENE
! this factor comes from $b_a$ argument in the Bessel. Kperp is not used otherwise.
kparray(iky, ikx, iz, eo) = &
- SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/ hatB(iz,eo)
+ SQRT( gxx(iz,eo)*kx**2 + 2._xp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/ hatB(iz,eo)
ENDDO
ENDDO
ENDDO
ENDDO
- two_third_kpmax = 2._dp/3._dp * MAXVAL(kparray)
+ two_third_kpmax = 2._xp/3._xp * MAXVAL(kparray)
END SUBROUTINE
SUBROUTINE grid_outputinputs(fid)
diff --git a/src/inital.F90 b/src/inital.F90
index d9813ba3..cfc88cb4 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -100,12 +100,12 @@ SUBROUTINE init_moments
ngp, ngj, ngz, iky0, contains_ky0, AA_x, AA_y
USE initial_par,ONLY: iseed, init_noiselvl, init_background
USE fields, ONLY: moments
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
USE model, ONLY: LINEARITY
USE parallel, ONLY: my_id
IMPLICIT NONE
- REAL(dp) :: noise
+ REAL(xp) :: noise
INTEGER, DIMENSION(12) :: iseedarr
INTEGER :: ia,ip,ij,ikx,iky,iz, ipi,iji,izi
@@ -124,7 +124,7 @@ SUBROUTINE init_moments
DO iz=1,local_nz
izi = iz+ngz/2
CALL RANDOM_NUMBER(noise)
- moments(ia,ipi,iji,iky,ikx,izi,:) = (init_background + init_noiselvl*(noise-0.5_dp))
+ moments(ia,ipi,iji,iky,ikx,izi,:) = (init_background + init_noiselvl*(noise-0.5_xp))
END DO
END DO
END DO
@@ -163,20 +163,20 @@ SUBROUTINE init_gyrodens
USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
ngp, ngj, ngz, iky0, parray, jarray, contains_ky0, AA_x, AA_y
USE fields, ONLY: moments
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
USE initial_par,ONLY: iseed, init_noiselvl, init_background
USE model, ONLY: LINEARITY
USE parallel, ONLY: my_id
IMPLICIT NONE
- REAL(dp) :: noise
+ REAL(xp) :: noise
INTEGER :: ia,ip,ij,ikx,iky,iz
INTEGER, DIMENSION(12) :: iseedarr
! Seed random number generator
iseedarr(:)=iseed
CALL RANDOM_SEED(PUT=iseedarr+my_id)
- moments = 0._dp
+ moments = 0._xp
!**** Broad noise initialization *******************************************
DO ia=1,local_na
DO ip=1+ngp/2,local_np+ngp/2
@@ -186,9 +186,9 @@ SUBROUTINE init_gyrodens
DO iz=1+ngz/2,local_nz+ngz/2
CALL RANDOM_NUMBER(noise)
IF ( (parray(ip) .EQ. 0) .AND. (jarray(ij) .EQ. 0) ) THEN
- moments(ia,ip,ij,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
+ moments(ia,ip,ij,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_xp))
ELSE
- moments(ia,ip,ij,iky,ikx,iz,:) = 0._dp
+ moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
ENDIF
END DO
END DO
@@ -225,12 +225,12 @@ SUBROUTINE init_phi
USE grid, ONLY: total_nkx, local_nky, local_nz,&
ngz, iky0, ikx0, contains_ky0
USE fields, ONLY: phi, moments
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
USE initial_par,ONLY: iseed, init_noiselvl, init_background
USE model, ONLY: LINEARITY
USE parallel, ONLY: my_id
IMPLICIT NONE
- REAL(dp) :: noise
+ REAL(xp) :: noise
INTEGER, DIMENSION(12) :: iseedarr
INTEGER :: iky,ikx,iz
! Seed random number generator
@@ -241,7 +241,7 @@ SUBROUTINE init_phi
DO iky=1,local_nky
DO iz=1,local_nz+ngz
CALL RANDOM_NUMBER(noise)
- phi(iky,ikx,iz) = (init_background + init_noiselvl*(noise-0.5_dp))!*AA_x(ikx)*AA_y(iky)
+ phi(iky,ikx,iz) = (init_background + init_noiselvl*(noise-0.5_xp))!*AA_x(ikx)*AA_y(iky)
ENDDO
END DO
END DO
@@ -251,11 +251,11 @@ SUBROUTINE init_phi
DO ikx=2,total_nkx/2
phi(iky0,ikx,iz) = phi(iky0,total_nkx+2-ikx,iz)
ENDDO
- phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz),dp) !origin must be real
+ phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz),xp) !origin must be real
END DO
ENDIF
!**** ensure no previous moments initialization
- moments = 0._dp
+ moments = 0._xp
!**** Zonal Flow initialization *******************************************
! put a mode at ikx = mode number + 1, symmetry is already included since kx>=0
! IF(INIT_ZF .GT. 0) THEN
@@ -263,8 +263,8 @@ SUBROUTINE init_phi
! IF( (INIT_ZF+1 .GT. ikxs) .AND. (INIT_ZF+1 .LT. ikxe) ) THEN
! DO ia=1,local_na
! DO iz = 1,local_nz+ngz
- ! phi(iky0,INIT_ZF+1,iz) = ZF_AMP*(2._dp*PI)**2/deltakx/deltaky/2._dp * COS((iz-1)/Nz*2._dp*PI)
- ! moments(ia,ip0,ij0,iky0,INIT_ZF+1,iz,:) = kxarray(INIT_ZF+1)**2*phi(iky0,INIT_ZF+1,iz)* COS((iz-1)/Nz*2._dp*PI)
+ ! phi(iky0,INIT_ZF+1,iz) = ZF_AMP*(2._xp*PI)**2/deltakx/deltaky/2._xp * COS((iz-1)/Nz*2._xp*PI)
+ ! moments(ia,ip0,ij0,iky0,INIT_ZF+1,iz,:) = kxarray(INIT_ZF+1)**2*phi(iky0,INIT_ZF+1,iz)* COS((iz-1)/Nz*2._xp*PI)
! ENDDO
! ENDDO
! ENDIF
@@ -279,14 +279,14 @@ SUBROUTINE init_phi_ppj
USE grid, ONLY: total_nkx, local_nky, local_nz,&
ngz, iky0, ikx0, contains_ky0, ieven, kxarray, kyarray, zarray, deltakx
USE fields, ONLY: phi, moments
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
USE initial_par,ONLY: iseed, init_noiselvl, init_background
USE model, ONLY: LINEARITY
USE geometry, ONLY: Jacobian, iInt_Jacobian
IMPLICIT NONE
- REAL(dp) :: kx, ky, z, amp
+ REAL(xp) :: kx, ky, z, amp
INTEGER :: ikx, iky, iz
- amp = 1.0_dp
+ amp = 1.0_xp
!**** ppj initialization *******************************************
DO ikx=1,total_nkx
kx = kxarray(ikx)
@@ -295,9 +295,9 @@ SUBROUTINE init_phi_ppj
DO iz=1,local_nz+ngz
z = zarray(iz,ieven)
IF (ky .NE. 0) THEN
- phi(iky,ikx,iz) = 0._dp
+ phi(iky,ikx,iz) = 0._xp
ELSE
- phi(iky,ikx,iz) = 0.5_dp*amp*(deltakx/(ABS(kx)+deltakx))
+ phi(iky,ikx,iz) = 0.5_xp*amp*(deltakx/(ABS(kx)+deltakx))
ENDIF
! z-dep and noise
phi(iky,ikx,iz) = phi(iky,ikx,iz) * &
@@ -311,11 +311,11 @@ SUBROUTINE init_phi_ppj
DO ikx=2,total_nkx/2
phi(iky0,ikx,iz) = phi(iky0,total_nkx+2-ikx,iz)
ENDDO
- phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz),dp) !origin must be real
+ phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz),xp) !origin must be real
END DO
ENDIF
!**** ensure no previous moments initialization
- moments = 0._dp
+ moments = 0._xp
END SUBROUTINE init_phi_ppj
!******************************************************************************!
@@ -328,12 +328,12 @@ SUBROUTINE initialize_blob
AA_x, AA_y, parray, jarray,&
ngp,ngj,ngz, iky0, ieven, kxarray, kyarray, zarray
USE fields, ONLY: moments
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
USE initial_par,ONLY: iseed, init_noiselvl, init_background
USE model, ONLY: LINEARITY
USE geometry, ONLY: Jacobian, iInt_Jacobian, shear
IMPLICIT NONE
- REAL(dp) ::kx, ky, z, sigma_x, sigma_y, gain
+ REAL(xp) ::kx, ky, z, sigma_x, sigma_y, gain
INTEGER :: ia,iky,ikx,iz,ip,ij, p, j
sigma_y = 1.0
sigma_x = sigma_y
@@ -372,15 +372,15 @@ SUBROUTINE init_ppj
AA_x, AA_y, deltakx, deltaky,contains_ky0,&
ngp,ngj,ngz, iky0, ieven, kxarray, kyarray, zarray
USE fields, ONLY: moments
- USE prec_const, ONLY: dp, pi
+ USE prec_const, ONLY: xp, pi
USE initial_par,ONLY: iseed, init_noiselvl, init_background
USE model, ONLY: LINEARITY
USE geometry, ONLY: Jacobian, iInt_Jacobian, shear
IMPLICIT NONE
- REAL(dp) :: kx, ky, sigma_z, amp, z
+ REAL(xp) :: kx, ky, sigma_z, amp, z
INTEGER :: ia,iky,ikx,iz,ip,ij
- sigma_z = pi/4._dp
- amp = 1.0_dp
+ sigma_z = pi/4._xp
+ amp = 1.0_xp
!**** Broad noise initialization *******************************************
DO ia=1,local_na
DO ip=1,local_np+ngp
@@ -394,15 +394,15 @@ SUBROUTINE init_ppj
z = zarray(iz,ieven)
IF (kx .EQ. 0) THEN
IF(ky .EQ. 0) THEN
- moments(ia,ip,ij,iky,ikx,iz,:) = 0._dp
+ moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
ELSE
- moments(ia,ip,ij,iky,ikx,iz,:) = 0.5_dp * deltaky/(ABS(ky)+deltaky)
+ moments(ia,ip,ij,iky,ikx,iz,:) = 0.5_xp * deltaky/(ABS(ky)+deltaky)
ENDIF
ELSE
IF(ky .GT. 0) THEN
moments(ia,ip,ij,iky,ikx,iz,:) = (deltakx/(ABS(kx)+deltakx))*(deltaky/(ABS(ky)+deltaky))
ELSE
- moments(ia,ip,ij,iky,ikx,iz,:) = 0.5_dp*amp*(deltakx/(ABS(kx)+deltakx))
+ moments(ia,ip,ij,iky,ikx,iz,:) = 0.5_xp*amp*(deltakx/(ABS(kx)+deltakx))
ENDIF
ENDIF
! z-dep and noise
@@ -419,7 +419,7 @@ SUBROUTINE init_ppj
END DO
ENDIF
ELSE
- moments(ia,ip,ij,:,:,:,:) = 0._dp
+ moments(ia,ip,ij,:,:,:,:) = 0._xp
ENDIF
END DO
END DO
diff --git a/src/initial_par_mod.F90 b/src/initial_par_mod.F90
index 0b95194f..4fd93b0b 100644
--- a/src/initial_par_mod.F90
+++ b/src/initial_par_mod.F90
@@ -9,13 +9,13 @@ MODULE initial_par
CHARACTER(len=32), PUBLIC, PROTECTED :: INIT_OPT = 'phi'
! Initialization through a zonal flow phi
INTEGER, PUBLIC, PROTECTED :: INIT_ZF = 0
- REAL(DP), PUBLIC, PROTECTED :: ZF_AMP = 1E+3_dp
+ REAL(xp), PUBLIC, PROTECTED :: ZF_AMP = 1E+3_xp
! Act on modes artificially (keep/wipe, zonal, non zonal, entropy mode etc.)
CHARACTER(len=32), PUBLIC, PROTECTED :: ACT_ON_MODES = 'nothing'
! Initial background level
- REAL(dp), PUBLIC, PROTECTED :: init_background=0._dp
+ REAL(xp), PUBLIC, PROTECTED :: init_background=0._xp
! Initial noise amplitude
- REAL(dp), PUBLIC, PROTECTED :: init_noiselvl=1E-6_dp
+ REAL(xp), PUBLIC, PROTECTED :: init_noiselvl=1E-6_xp
! Initialization for random number generator
INTEGER, PUBLIC, PROTECTED :: iseed=42
diff --git a/src/lag_interp_mod.F90 b/src/lag_interp_mod.F90
index 8415df5d..39cea31f 100644
--- a/src/lag_interp_mod.F90
+++ b/src/lag_interp_mod.F90
@@ -21,11 +21,11 @@ CONTAINS
!> Third order lagrange interpolation
SUBROUTINE lag3interp_scalar(y_in,x_in,n_in,y_out,x_out)
INTEGER, INTENT(IN) :: n_in
- REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
- REAL(dp), INTENT(IN) :: x_out
- REAL(dp), INTENT(OUT) :: y_out
+ REAL(xp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
+ REAL(xp), INTENT(IN) :: x_out
+ REAL(xp), INTENT(OUT) :: y_out
- REAL(dp), DIMENSION(1) :: xout_wrap, yout_wrap
+ REAL(xp), DIMENSION(1) :: xout_wrap, yout_wrap
xout_wrap = x_out
call lag3interp_array(y_in,x_in,n_in,yout_wrap,xout_wrap,1)
@@ -36,11 +36,11 @@ CONTAINS
!> Third order lagrange interpolation
subroutine lag3interp_array(y_in,x_in,n_in,y_out,x_out,n_out)
INTEGER, INTENT(IN) :: n_in,n_out
- REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
- REAL(dp), DIMENSION(n_out), INTENT(IN) :: x_out
- REAL(dp), DIMENSION(n_out), INTENT(OUT) :: y_out
+ REAL(xp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
+ REAL(xp), DIMENSION(n_out), INTENT(IN) :: x_out
+ REAL(xp), DIMENSION(n_out), INTENT(OUT) :: y_out
- REAL(dp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
+ REAL(xp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
INTEGER :: j,jm,j0,j1,j2
INTEGER :: jstart,jfirst,jlast,jstep
@@ -91,11 +91,11 @@ CONTAINS
SUBROUTINE lag3interp_complex(y_in,x_in,n_in,y_out,x_out,n_out)
INTEGER, INTENT(IN) :: n_in,n_out
COMPLEX, DIMENSION(n_in), INTENT(IN) :: y_in
- REAL(dp), DIMENSION(n_in), INTENT(IN) :: x_in
+ REAL(xp), DIMENSION(n_in), INTENT(IN) :: x_in
COMPLEX, DIMENSION(n_out), INTENT(OUT) :: y_out
- REAL(dp), DIMENSION(n_out), INTENT(IN) :: x_out
+ REAL(xp), DIMENSION(n_out), INTENT(IN) :: x_out
- REAL(dp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
+ REAL(xp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
INTEGER :: j,jm,j0,j1,j2
INTEGER :: jstart,jfirst,jlast,jstep
@@ -143,22 +143,22 @@ CONTAINS
!>2D interpolation
- !\TODO check whether a "REAL(dp)" 2D interpolation would
+ !\TODO check whether a "REAL(xp)" 2D interpolation would
!! be more appropriate
SUBROUTINE lag3interp_2d(y_in,x1_in,n1_in,x2_in,n2_in,&
&y_out,x1_out,n1_out,x2_out,n2_out)
INTEGER, INTENT(IN) :: n1_in,n2_in,n1_out,n2_out
- REAL(dp), DIMENSION(n1_in,n2_in), INTENT(IN) :: y_in
- REAL(dp), DIMENSION(n1_in) :: x1_in
- REAL(dp), DIMENSION(n2_in) :: x2_in
- REAL(dp), DIMENSION(n1_out), INTENT(IN) :: x1_out
- REAL(dp), DIMENSION(n2_out), INTENT(IN) :: x2_out
- REAL(dp), DIMENSION(n1_out,n2_out), INTENT(OUT) :: y_out
+ REAL(xp), DIMENSION(n1_in,n2_in), INTENT(IN) :: y_in
+ REAL(xp), DIMENSION(n1_in) :: x1_in
+ REAL(xp), DIMENSION(n2_in) :: x2_in
+ REAL(xp), DIMENSION(n1_out), INTENT(IN) :: x1_out
+ REAL(xp), DIMENSION(n2_out), INTENT(IN) :: x2_out
+ REAL(xp), DIMENSION(n1_out,n2_out), INTENT(OUT) :: y_out
!local variables
- REAL(dp), DIMENSION(n2_in) :: y2_in_tmp
- REAL(dp), DIMENSION(n2_out) :: y2_out_tmp
- REAL(dp), DIMENSION(n1_in,n2_out) :: y_tmp
+ REAL(xp), DIMENSION(n2_in) :: y2_in_tmp
+ REAL(xp), DIMENSION(n2_out) :: y2_out_tmp
+ REAL(xp), DIMENSION(n1_in,n2_out) :: y_tmp
INTEGER :: i
DO i=1,n1_in
@@ -181,11 +181,11 @@ CONTAINS
IMPLICIT NONE
INTEGER, INTENT(IN) :: n_in
- REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
- REAL(dp), INTENT(IN) :: x_out
- REAL(dp), INTENT(OUT) :: dydx_out
+ REAL(xp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
+ REAL(xp), INTENT(IN) :: x_out
+ REAL(xp), INTENT(OUT) :: dydx_out
- REAL(dp), DIMENSION(1) :: xout_wrap, dydxout_wrap
+ REAL(xp), DIMENSION(1) :: xout_wrap, dydxout_wrap
xout_wrap = x_out
call lag3deriv_array(y_in,x_in,n_in,dydxout_wrap,xout_wrap,1)
@@ -198,11 +198,11 @@ CONTAINS
!>Returns Derivative based on a 3rd order lagrange interpolation
subroutine lag3deriv_array(y_in,x_in,n_in,dydx_out,x_out,n_out)
INTEGER :: n_in,n_out
- REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
- REAL(dp), DIMENSION(n_out), INTENT(IN) :: x_out
- REAL(dp), DIMENSION(n_out), INTENT(OUT) :: dydx_out
+ REAL(xp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
+ REAL(xp), DIMENSION(n_out), INTENT(IN) :: x_out
+ REAL(xp), DIMENSION(n_out), INTENT(OUT) :: dydx_out
- REAL(dp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
+ REAL(xp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
INTEGER :: j,jm,j0,j1,j2
INTEGER :: jstart,jfirst,jlast,jstep
diff --git a/src/miller_mod.F90 b/src/miller_mod.F90
index 3ecce924..495347eb 100644
--- a/src/miller_mod.F90
+++ b/src/miller_mod.F90
@@ -20,15 +20,15 @@ MODULE miller
private
- real(dp) :: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta
- real(dp) :: thetaShift
- real(dp) :: thetak, thetad
+ real(xp) :: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta
+ real(xp) :: thetaShift
+ real(xp) :: thetak, thetad
INTEGER :: ierr
CONTAINS
!>Set defaults for miller parameters
subroutine set_miller_parameters(kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_)
- real(dp), INTENT(IN) :: kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_
+ real(xp), INTENT(IN) :: kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_
rho = -1.0
kappa = kappa_
s_kappa = s_kappa_
@@ -46,62 +46,62 @@ CONTAINS
!>Get Miller metric, magnetic field, jacobian etc.
subroutine get_miller(trpeps,major_R,major_Z,q0,shat,Npol,amhd,edge_opt,&
- C_y,C_xy,dpdx_pm_geom,gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,dBdx_,dBdy_,&
+ C_y,C_xy,xpdx_pm_geom,gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,dBdx_,dBdy_,&
Bfield_,jacobian_,dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_,Ckxky_,gradz_coeff_)
!!!!!!!!!!!!!!!! GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real(dp), INTENT(INOUT) :: trpeps ! eps in gyacomo (inverse aspect ratio)
- real(dp), INTENT(INOUT) :: major_R ! major radius
- real(dp), INTENT(INOUT) :: major_Z ! major Z
- real(dp), INTENT(INOUT) :: q0 ! safetyfactor
- real(dp), INTENT(INOUT) :: shat ! safetyfactor
+ real(xp), INTENT(INOUT) :: trpeps ! eps in gyacomo (inverse aspect ratio)
+ real(xp), INTENT(INOUT) :: major_R ! major radius
+ real(xp), INTENT(INOUT) :: major_Z ! major Z
+ real(xp), INTENT(INOUT) :: q0 ! safetyfactor
+ real(xp), INTENT(INOUT) :: shat ! safetyfactor
INTEGER, INTENT(IN) :: Npol ! number of poloidal turns
- real(dp), INTENT(INOUT) :: amhd ! alpha mhd
- real(dp), INTENT(INOUT) :: edge_opt ! alpha mhd
- real(dp), INTENT(INOUT) :: dpdx_pm_geom ! amplitude mag. eq. pressure grad.
- real(dp), INTENT(INOUT) :: C_y, C_xy
- real(dp), dimension(1:local_nz+Ngz,1:2), INTENT(INOUT) :: &
+ real(xp), INTENT(INOUT) :: amhd ! alpha mhd
+ real(xp), INTENT(INOUT) :: edge_opt ! alpha mhd
+ real(xp), INTENT(INOUT) :: xpdx_pm_geom ! amplitude mag. eq. pressure grad.
+ real(xp), INTENT(INOUT) :: C_y, C_xy
+ real(xp), dimension(1:local_nz+Ngz,1:2), INTENT(INOUT) :: &
gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,&
dBdx_,dBdy_,Bfield_,jacobian_,&
dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_, &
gradz_coeff_
- real(dp), dimension(1:local_Nky,1:local_Nkx,1:local_nz+Ngz,1:2), INTENT(INOUT) :: Ckxky_
+ real(xp), dimension(1:local_Nky,1:local_Nkx,1:local_nz+Ngz,1:2), INTENT(INOUT) :: Ckxky_
INTEGER :: iz, ikx, iky, eo
! No parameter in gyacomo yet
- real(dp) :: sign_Ip_CW=1 ! current sign (only normal current)
- real(dp) :: sign_Bt_CW=1 ! current sign (only normal current)
+ real(xp) :: sign_Ip_CW=1 ! current sign (only normal current)
+ real(xp) :: sign_Bt_CW=1 ! current sign (only normal current)
!!!!!!!!!!!!!! END GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Auxiliary variables for curvature computation
- real(dp) :: G1,G2,G3,Cx,Cy,ky,kx
+ real(xp) :: G1,G2,G3,Cx,Cy,ky,kx
integer:: np, np_s, Npol_ext, Npol_s
- real(dp), dimension(500*(Npol+2)):: R,Z,R_rho,Z_rho,R_theta,Z_theta,R_theta_theta,Z_theta_theta,dlp,Rc,cosu,sinu,Bphi
- real(dp), dimension(500*(Npol+2)):: drRcirc, drRelong, drRelongTilt, drRtri, drRtriTilt, drZcirc, drZelong, drZelongTilt
- real(dp), dimension(500*(Npol+2)):: drZtri, drZtriTilt, dtdtRcirc, dtdtRelong, dtdtRelongTilt, dtdtRtri, dtdtRtriTilt
- real(dp), dimension(500*(Npol+2)):: dtdtZcirc, dtdtZelong, dtdtZelongTilt, dtdtZtri, dtdtZtriTilt, dtRcirc, dtRelong
- real(dp), dimension(500*(Npol+2)):: dtRelongTilt, dtRtri, dtRtriTilt, dtZcirc, dtZelong, dtZelongTilt, dtZtri, dtZtriTilt
- real(dp), dimension(500*(Npol+2)):: Rcirc, Relong, RelongTilt, Rtri, RtriTilt, Zcirc, Zelong, ZelongTilt, Ztri, ZtriTilt
- ! real(dp), dimension(500*(Npol+2)):: drrShape, drrAng, drxAng, dryAng, dtdtrShape, dtdtrAng, dtdtxAng
- ! real(dp), dimension(500*(Npol+2)):: dtdtyAng, dtrShape, dtrAng, dtxAng, dtyAng, rShape, rAng, xAng, yAng
- real(dp), dimension(500*(Npol+2)):: theta, thAdj, J_r, B, Bp, D0, D1, D2, D3, nu, chi, psi1, nu1
- real(dp), dimension(500*(Npol+2)):: tmp_reverse, theta_reverse, tmp_arr
-
- real(dp), dimension(500*(Npol+1)):: theta_s, thAdj_s, chi_s, theta_s_reverse
- real(dp), dimension(500*(Npol+1)):: R_s, Z_s, R_theta_s, Z_theta_s, Rc_s, cosu_s, sinu_s, Bphi_s, B_s, Bp_s, dlp_s
- real(dp), dimension(500*(Npol+1)):: dtRcirc_s, dtRelong_s, dtRelongTilt_s, dtRtri_s, dtRtriTilt_s, dtZcirc_s
- real(dp), dimension(500*(Npol+1)):: dtZelong_s, dtZelongTilt_s, dtZtri_s, dtZtriTilt_s, Rcirc_s, Relong_s, RelongTilt_s
- real(dp), dimension(500*(Npol+1)):: Rtri_s, RtriTilt_s, Zcirc_s, Zelong_s, ZelongTilt_s, Ztri_s, ZtriTilt_s!, dtrShape_s
- ! real(dp), dimension(500*(Npol+1)):: dtrAng_s, dtxAng_s, dtyAng_s, rShape_s, rAng_s, xAng_s, yAng_s
- real(dp), dimension(500*(Npol+1)):: psi1_s, nu1_s, dchidx_s, dB_drho_s, dB_dl_s, dnu_drho_s, dnu_dl_s, grad_nu_s
- real(dp), dimension(500*(Npol+1)):: gxx, gxy, gxz, gyy, gyz, gzz, dtheta_dchi_s, dBp_dchi_s, jacobian, dBdx, dBdz
- real(dp), dimension(500*(Npol+1)):: g_xx, g_xy, g_xz, g_yy, g_yz, g_zz, tmp_arr_s, dxdR_s, dxdZ_s, K_x, K_y !tmp_arr2
-
- real(dp), dimension(1:Nz):: gxx_out,gxy_out,gxz_out,gyy_out,gyz_out,gzz_out,Bfield_out,jacobian_out, dBdx_out, dBdz_out, chi_out
- real(dp), dimension(1:Nz):: R_out, Z_out, dxdR_out, dxdZ_out
- real(dp):: d_inv, drPsi, dxPsi, dq_dx, dq_dpsi, R0, Z0, B0, F, D0_full, D1_full, D2_full, D3_full
- !real(dp) :: Lnorm, Psi0 ! currently module-wide defined anyway
- real(dp):: pprime, ffprime, D0_mid, D1_mid, D2_mid, D3_mid, dx_drho, pi, mu_0, dzprimedz
- ! real(dp):: rho_a, psiN, drpsiN, CN2, CN3, Rcenter, Zcenter, drRcenter, drZcenter
+ real(xp), dimension(500*(Npol+2)):: R,Z,R_rho,Z_rho,R_theta,Z_theta,R_theta_theta,Z_theta_theta,dlp,Rc,cosu,sinu,Bphi
+ real(xp), dimension(500*(Npol+2)):: drRcirc, drRelong, drRelongTilt, drRtri, drRtriTilt, drZcirc, drZelong, drZelongTilt
+ real(xp), dimension(500*(Npol+2)):: drZtri, drZtriTilt, dtdtRcirc, dtdtRelong, dtdtRelongTilt, dtdtRtri, dtdtRtriTilt
+ real(xp), dimension(500*(Npol+2)):: dtdtZcirc, dtdtZelong, dtdtZelongTilt, dtdtZtri, dtdtZtriTilt, dtRcirc, dtRelong
+ real(xp), dimension(500*(Npol+2)):: dtRelongTilt, dtRtri, dtRtriTilt, dtZcirc, dtZelong, dtZelongTilt, dtZtri, dtZtriTilt
+ real(xp), dimension(500*(Npol+2)):: Rcirc, Relong, RelongTilt, Rtri, RtriTilt, Zcirc, Zelong, ZelongTilt, Ztri, ZtriTilt
+ ! real(xp), dimension(500*(Npol+2)):: drrShape, drrAng, drxAng, dryAng, dtdtrShape, dtdtrAng, dtdtxAng
+ ! real(xp), dimension(500*(Npol+2)):: dtdtyAng, dtrShape, dtrAng, dtxAng, dtyAng, rShape, rAng, xAng, yAng
+ real(xp), dimension(500*(Npol+2)):: theta, thAdj, J_r, B, Bp, D0, D1, D2, D3, nu, chi, psi1, nu1
+ real(xp), dimension(500*(Npol+2)):: tmp_reverse, theta_reverse, tmp_arr
+
+ real(xp), dimension(500*(Npol+1)):: theta_s, thAdj_s, chi_s, theta_s_reverse
+ real(xp), dimension(500*(Npol+1)):: R_s, Z_s, R_theta_s, Z_theta_s, Rc_s, cosu_s, sinu_s, Bphi_s, B_s, Bp_s, dlp_s
+ real(xp), dimension(500*(Npol+1)):: dtRcirc_s, dtRelong_s, dtRelongTilt_s, dtRtri_s, dtRtriTilt_s, dtZcirc_s
+ real(xp), dimension(500*(Npol+1)):: dtZelong_s, dtZelongTilt_s, dtZtri_s, dtZtriTilt_s, Rcirc_s, Relong_s, RelongTilt_s
+ real(xp), dimension(500*(Npol+1)):: Rtri_s, RtriTilt_s, Zcirc_s, Zelong_s, ZelongTilt_s, Ztri_s, ZtriTilt_s!, dtrShape_s
+ ! real(xp), dimension(500*(Npol+1)):: dtrAng_s, dtxAng_s, dtyAng_s, rShape_s, rAng_s, xAng_s, yAng_s
+ real(xp), dimension(500*(Npol+1)):: psi1_s, nu1_s, dchidx_s, dB_drho_s, dB_dl_s, dnu_drho_s, dnu_dl_s, grad_nu_s
+ real(xp), dimension(500*(Npol+1)):: gxx, gxy, gxz, gyy, gyz, gzz, dtheta_dchi_s, dBp_dchi_s, jacobian, dBdx, dBdz
+ real(xp), dimension(500*(Npol+1)):: g_xx, g_xy, g_xz, g_yy, g_yz, g_zz, tmp_arr_s, dxdR_s, dxdZ_s, K_x, K_y !tmp_arr2
+
+ real(xp), dimension(1:Nz):: gxx_out,gxy_out,gxz_out,gyy_out,gyz_out,gzz_out,Bfield_out,jacobian_out, dBdx_out, dBdz_out, chi_out
+ real(xp), dimension(1:Nz):: R_out, Z_out, dxdR_out, dxdZ_out
+ real(xp):: d_inv, drPsi, dxPsi, dq_dx, dq_xpsi, R0, Z0, B0, F, D0_full, D1_full, D2_full, D3_full
+ !real(xp) :: Lnorm, Psi0 ! currently module-wide defined anyway
+ real(xp):: pprime, ffprime, D0_mid, D1_mid, D2_mid, D3_mid, dx_drho, pi, mu_0, dzprimedz
+ ! real(xp):: rho_a, psiN, drpsiN, CN2, CN3, Rcenter, Zcenter, drRcenter, drZcenter
logical:: bMaxShift
integer:: i, k, iBmax
@@ -123,7 +123,7 @@ CONTAINS
pi = acos(-1.0)
mu_0=4.0*pi
- theta=linspace(-pi*Npol_ext,pi*Npol_ext-2._dp*pi*Npol_ext/np,np)
+ theta=linspace(-pi*Npol_ext,pi*Npol_ext-2._xp*pi*Npol_ext/np,np)
d_inv=asin(delta)
thetaShift = 0.0
@@ -263,11 +263,11 @@ CONTAINS
!--------shear is expected to be defined as rho/q*dq/drho--------!
dq_dx=shat*q0/rho/dx_drho
- dq_dpsi=dq_dx/dxPsi
+ dq_xpsi=dq_dx/dxPsi
pprime=-amhd/q0**2/R0/(2*mu_0)*B0**2/drPsi
- !neg. dpdx normalized to magnetic pressure for pressure term
- dpdx_pm_geom=amhd/q0**2/R0/dx_drho
+ !neg. xpdx normalized to magnetic pressure for pressure term
+ xpdx_pm_geom=amhd/q0**2/R0/dx_drho
!first coefficient of psi in varrho expansion
psi1 = R*Bp*sign_Ip_CW
@@ -296,7 +296,7 @@ CONTAINS
D2_mid=D2(np/2+1)
D3_mid=D3(np/2+1)
- ffprime=-(sign_Ip_CW*dq_dpsi*2.*pi*Npol_ext+D0_full+D2_full*pprime)/D1_full
+ ffprime=-(sign_Ip_CW*dq_xpsi*2.*pi*Npol_ext+D0_full+D2_full*pprime)/D1_full
if (my_id==0) then
write(*,'(A,ES12.4)') "ffprime = ", ffprime
@@ -405,10 +405,10 @@ CONTAINS
!contravariant metric coefficients (varrho,l,fz)->(x,y,z)
gxx=(psi1_s/dxPsi)**2
gxy=-psi1_s/dxPsi*C_y*sign_Ip_CW*nu1_s
- gxz=-psi1_s/dxPsi*(nu1_s+psi1_s*dq_dpsi*chi_s)/q0
+ gxz=-psi1_s/dxPsi*(nu1_s+psi1_s*dq_xpsi*chi_s)/q0
gyy=C_y**2*(grad_nu_s**2+1/R_s**2)
- gyz=sign_Ip_CW*C_y/q0*(grad_nu_s**2+dq_dpsi*nu1_s*psi1_s*chi_s)
- gzz=1./q0**2*(grad_nu_s**2+2.*dq_dpsi*nu1_s*psi1_s*chi_s+(dq_dpsi*psi1_s*chi_s)**2)
+ gyz=sign_Ip_CW*C_y/q0*(grad_nu_s**2+dq_xpsi*nu1_s*psi1_s*chi_s)
+ gzz=1./q0**2*(grad_nu_s**2+2.*dq_xpsi*nu1_s*psi1_s*chi_s+(dq_xpsi*psi1_s*chi_s)**2)
jacobian=1./sqrt(gxx*gyy*gzz + 2.*gxy*gyz*gxz - gxz**2*gyy - gyz**2*gxx - gzz*gxy**2)
@@ -422,7 +422,7 @@ CONTAINS
!Bfield derivatives
!dBdx = e_x * nabla B = J (nabla y x nabla z) * nabla B
- dBdx=jacobian*C_y/(q0*R_s)*(F/(R_s*psi1_s)*dB_drho_s+(nu1_s+dq_dpsi*chi_s*psi1_s)*dB_dl_s)
+ dBdx=jacobian*C_y/(q0*R_s)*(F/(R_s*psi1_s)*dB_drho_s+(nu1_s+dq_xpsi*chi_s*psi1_s)*dB_dl_s)
dBdz=1./B_s*(Bp_s*dBp_dchi_s-F**2/R_s**3*R_theta_s)
!curvature terms (these are just local and will be recalculated in geometry.F90)
@@ -526,7 +526,7 @@ CONTAINS
ENDDO
ENDDO
! coefficient in the front of parallel derivative
- gradz_coeff_(iz,eo) = 1._dp / jacobian_(iz,eo) / Bfield_(iz,eo)
+ gradz_coeff_(iz,eo) = 1._xp / jacobian_(iz,eo) / Bfield_(iz,eo)
ENDDO
ENDDO
@@ -536,8 +536,8 @@ CONTAINS
SUBROUTINE update_ghosts_z(fz_)
IMPLICIT NONE
! INTEGER, INTENT(IN) :: nztot_
- REAL(dp), DIMENSION(1:local_nz+Ngz), INTENT(INOUT) :: fz_
- REAL(dp), DIMENSION(-2:2) :: buff
+ REAL(xp), DIMENSION(1:local_nz+Ngz), INTENT(INOUT) :: fz_
+ REAL(xp), DIMENSION(-2:2) :: buff
INTEGER :: status(MPI_STATUS_SIZE), count, last, first
last = local_nz+Ngz/2
first= 1 + Ngz/2
@@ -578,9 +578,9 @@ CONTAINS
!> Generate an equidistant array from min to max with n points
function linspace(min,max,n) result(out)
- real(dp), INTENT(IN):: min, max
+ real(xp), INTENT(IN):: min, max
integer, INTENT(IN):: n
- real(dp), dimension(n):: out
+ real(xp), dimension(n):: out
do i=1,n
out(i)=min+(i-1)*(max-min)/(n-1)
@@ -588,20 +588,20 @@ CONTAINS
end function linspace
!> Weighted average
- real(dp) function average(var,weight)
- real(dp), dimension(np), INTENT(IN):: var, weight
+ real(xp) function average(var,weight)
+ real(xp), dimension(np), INTENT(IN):: var, weight
average=sum(var*weight)/sum(weight)
end function average
!> full theta integral with weight function dlp
- real(dp) function dlp_int(var,dlp)
- real(dp), dimension(np), INTENT(IN):: var, dlp
+ real(xp) function dlp_int(var,dlp)
+ real(xp), dimension(np), INTENT(IN):: var, dlp
dlp_int=sum(var*dlp)*2*pi*Npol_ext/np
end function dlp_int
!> theta integral with weight function dlp, up to index 'ind'
- real(dp) function dlp_int_ind(var,dlp,ind)
- real(dp), dimension(np), INTENT(IN):: var, dlp
+ real(xp) function dlp_int_ind(var,dlp,ind)
+ real(xp), dimension(np), INTENT(IN):: var, dlp
integer, INTENT(IN):: ind
dlp_int_ind=0.
@@ -615,8 +615,8 @@ CONTAINS
!> 1st derivative with 2nd order finite differences
function deriv_fd(y,x,n) result(out)
integer, INTENT(IN) :: n
- real(dp), dimension(n), INTENT(IN):: x,y
- real(dp), dimension(n) :: out,dx
+ real(xp), dimension(n), INTENT(IN):: x,y
+ real(xp), dimension(n) :: out,dx
!call lag3deriv(y,x,n,out,x,n)
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index ae9df467..6d93410a 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -9,23 +9,23 @@ MODULE model
INTEGER, PUBLIC, PROTECTED :: KERN = 0 ! Kernel model
CHARACTER(len=16), &
PUBLIC, PROTECTED ::LINEARITY= 'linear' ! To turn on non linear bracket term
- REAL(dp), PUBLIC, PROTECTED :: mu_x = 0._dp ! spatial x-Hyperdiffusivity coefficient (for num. stability)
- REAL(dp), PUBLIC, PROTECTED :: mu_y = 0._dp ! spatial y-Hyperdiffusivity coefficient (for num. stability)
+ REAL(xp), PUBLIC, PROTECTED :: mu_x = 0._xp ! spatial x-Hyperdiffusivity coefficient (for num. stability)
+ REAL(xp), PUBLIC, PROTECTED :: mu_y = 0._xp ! spatial y-Hyperdiffusivity coefficient (for num. stability)
INTEGER, PUBLIC, PROTECTED :: N_HD = 4 ! order of numerical spatial diffusion
LOGICAL, PUBLIC, PROTECTED :: HDz_h = .false. ! to apply z-hyperdiffusion on non adiab part
- REAL(dp), PUBLIC, PROTECTED :: mu_z = 0._dp ! spatial z-Hyperdiffusivity coefficient (for num. stability)
- REAL(dp), PUBLIC, PROTECTED :: mu_p = 0._dp ! kinetic para hyperdiffusivity coefficient (for num. stability)
- REAL(dp), PUBLIC, PROTECTED :: mu_j = 0._dp ! kinetic perp hyperdiffusivity coefficient (for num. stability)
+ REAL(xp), PUBLIC, PROTECTED :: mu_z = 0._xp ! spatial z-Hyperdiffusivity coefficient (for num. stability)
+ REAL(xp), PUBLIC, PROTECTED :: mu_p = 0._xp ! kinetic para hyperdiffusivity coefficient (for num. stability)
+ REAL(xp), PUBLIC, PROTECTED :: mu_j = 0._xp ! kinetic perp hyperdiffusivity coefficient (for num. stability)
CHARACTER(len=16), &
PUBLIC, PROTECTED :: HYP_V = 'hypcoll' ! hyperdiffusion model for velocity space ('none','hypcoll','dvpar4')
INTEGER, PUBLIC, PROTECTED :: Na = 1 ! number of evolved species
- REAL(dp), PUBLIC, PROTECTED :: nu = 0._dp ! collision frequency parameter
- REAL(dp), PUBLIC, PROTECTED :: k_gB = 1._dp ! Magnetic gradient strength (L_ref/L_gB)
- REAL(dp), PUBLIC, PROTECTED :: k_cB = 1._dp ! Magnetic curvature strength (L_ref/L_cB)
- REAL(dp), PUBLIC, PROTECTED :: lambdaD = 0._dp ! Debye length
- REAL(dp), PUBLIC, PROTECTED :: beta = 0._dp ! electron plasma Beta (8piNT_e/B0^2)
+ REAL(xp), PUBLIC, PROTECTED :: nu = 0._xp ! collision frequency parameter
+ REAL(xp), PUBLIC, PROTECTED :: k_gB = 1._xp ! Magnetic gradient strength (L_ref/L_gB)
+ REAL(xp), PUBLIC, PROTECTED :: k_cB = 1._xp ! Magnetic curvature strength (L_ref/L_cB)
+ REAL(xp), PUBLIC, PROTECTED :: lambdaD = 0._xp ! Debye length
+ REAL(xp), PUBLIC, PROTECTED :: beta = 0._xp ! electron plasma Beta (8piNT_e/B0^2)
LOGICAL, PUBLIC :: ADIAB_E = .false. ! adiabatic electron model
- REAL(dp), PUBLIC, PROTECTED :: tau_e = 1.0 ! electron temperature ratio for adiabatic electrons
+ REAL(xp), PUBLIC, PROTECTED :: tau_e = 1.0 ! electron temperature ratio for adiabatic electrons
! Auxiliary variable
LOGICAL, PUBLIC, PROTECTED :: EM = .false. ! Electromagnetic effects flag
PUBLIC :: model_readinputs, model_outputinputs
@@ -51,7 +51,7 @@ CONTAINS
IF(Na .EQ. 1) THEN
IF(my_id.EQ.0) print*, 'Adiabatic electron model -> beta = 0'
- beta = 0._dp
+ beta = 0._xp
ENDIF
IF(beta .GT. 0) THEN
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 81506404..c6efa288 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -15,22 +15,22 @@ SUBROUTINE compute_moments_eq_rhs
USE prec_const
USE collision
USE time_integration
- ! USE species, ONLY: dpdx
+ ! USE species, ONLY: xpdx
USE geometry, ONLY: gradz_coeff, dlnBdz, Ckxky!, Gamma_phipar
USE calculus, ONLY: interp_z, grad_z, grad_z2
- ! USE species, ONLY: dpdx
+ ! USE species, ONLY: xpdx
IMPLICIT NONE
INTEGER :: ia, iz, iky, ikx, ip ,ij, eo ! counters
INTEGER :: izi,ipi,iji ! interior points counters
INTEGER :: p_int, j_int ! loops indices and polynom. degrees
- REAL(dp) :: kx, ky, kperp2
- COMPLEX(dp) :: Tnapj, Tnapp2j, Tnapm2j, Tnapjp1, Tnapjm1 ! Terms from b x gradB and drives
- COMPLEX(dp) :: Tnapp1j, Tnapm1j, Tnapp1jm1, Tnapm1jm1 ! Terms from mirror force with non adiab moments_
- COMPLEX(dp) :: Ldamp, Fmir
- COMPLEX(dp) :: Mperp, Mpara, Dphi, Dpsi
- COMPLEX(dp) :: Unapm1j, Unapm1jp1, Unapm1jm1 ! Terms from mirror force with adiab moments_
- COMPLEX(dp) :: i_kx,i_ky
- COMPLEX(dp) :: Napj, RHS, phikykxz, psikykxz
+ REAL(xp) :: kx, ky, kperp2
+ COMPLEX(xp) :: Tnapj, Tnapp2j, Tnapm2j, Tnapjp1, Tnapjm1 ! Terms from b x gradB and drives
+ COMPLEX(xp) :: Tnapp1j, Tnapm1j, Tnapp1jm1, Tnapm1jm1 ! Terms from mirror force with non adiab moments_
+ COMPLEX(xp) :: Ldamp, Fmir
+ COMPLEX(xp) :: Mperp, Mpara, Dphi, Dpsi
+ COMPLEX(xp) :: Unapm1j, Unapm1jp1, Unapm1jm1 ! Terms from mirror force with adiab moments_
+ COMPLEX(xp) :: i_kx,i_ky
+ COMPLEX(xp) :: Napj, RHS, phikykxz, psikykxz
! Spatial loops
z:DO iz = 1,local_nz
izi = iz + ngz/2
@@ -54,7 +54,7 @@ SUBROUTINE compute_moments_eq_rhs
! Species loop
a:DO ia = 1,local_na
Napj = moments(ia,ipi,iji,iky,ikx,izi,updatetlevel)
- RHS = 0._dp
+ RHS = 0._xp
IF((CLOS .NE. 1) .OR. (p_int +2*j_int .LE. dmax)) THEN ! for the closure scheme
!! Compute moments_ mixing terms
! Perpendicular dynamic
@@ -95,7 +95,7 @@ SUBROUTINE compute_moments_eq_rhs
+xphijp1(ia,ip,ij)*kernel(ia,iji+1,iky,ikx,izi,eo) &
+xphijm1(ia,ip,ij)*kernel(ia,iji-1,iky,ikx,izi,eo) )*phikykxz
ELSE
- Dphi = 0._dp
+ Dphi = 0._xp
ENDIF
!! Vector potential term
IF ( (p_int .LE. 3) .AND. (p_int .GE. 1) ) THEN ! Kronecker p1 or p3
@@ -103,7 +103,7 @@ SUBROUTINE compute_moments_eq_rhs
+xpsijp1(ia,ip,ij)*kernel(ia,iji+1,iky,ikx,izi,eo) &
+xpsijm1(ia,ip,ij)*kernel(ia,iji-1,iky,ikx,izi,eo))*psikykxz
ELSE
- Dpsi = 0._dp
+ Dpsi = 0._xp
ENDIF
!! Sum of all RHS terms
RHS = &
@@ -118,7 +118,7 @@ SUBROUTINE compute_moments_eq_rhs
! Collision term
+ Capj(ia,ip,ij,iky,ikx,iz) &
! Perpendicular pressure effects (electromagnetic term) (TO CHECK)
- ! - i_ky*beta*dpdx(ia) * (Tnapj + Tnapp2j + Tnapm2j + Tnapjp1 + Tnapjm1)&
+ ! - i_ky*beta*xpdx(ia) * (Tnapj + Tnapp2j + Tnapm2j + Tnapjp1 + Tnapjm1)&
! Parallel drive term (should be negligible, to test)
! !! -Gamma_phipar(iz,eo)*Tphi*ddz_phi(iky,ikx,iz) &
! Numerical perpendicular hyperdiffusion
@@ -145,7 +145,7 @@ SUBROUTINE compute_moments_eq_rhs
CASE DEFAULT
END SELECT
ELSE
- RHS = 0._dp
+ RHS = 0._xp
ENDIF
!! Put RHS in the array
moments_rhs(ia,ip,ij,iky,ikx,iz,updatetlevel) = RHS
@@ -188,7 +188,7 @@ END SUBROUTINE compute_moments_eq_rhs
! USE grid, ONLY: contains_kx0, contains_ky0, ikx0, iky0,&
! ia,ias,iae,ip,ips,ipe, ij,ijs,ije, zarray,izs,ize
! IMPLICIT NONE
-! real(dp), DIMENSION(izs:ize) :: sinz
+! real(xp), DIMENSION(izs:ize) :: sinz
!
! sinz(izs:ize) = SIN(zarray(izs:ize,0))
!
@@ -201,19 +201,19 @@ END SUBROUTINE compute_moments_eq_rhs
! SELECT CASE (ip-1)
! CASE(0) ! Na00 term
! moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel) = moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel)&
-! +tau_q(ia) * sinz(izs:ize) * (1.5_dp*k_N(ia) - 1.125_dp*k_T(ia))
+! +tau_q(ia) * sinz(izs:ize) * (1.5_xp*k_N(ia) - 1.125_xp*k_T(ia))
! CASE(2) ! Na20 term
! moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel) = moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel)&
-! +tau_q(ia) * sinz(izs:ize) * (SQRT2*0.5_dp*k_N(ia) - 2.75_dp*k_T(ia))
+! +tau_q(ia) * sinz(izs:ize) * (SQRT2*0.5_xp*k_N(ia) - 2.75_xp*k_T(ia))
! CASE(4) ! Na40 term
! moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel) = moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel)&
-! +tau_q(ia) * sinz(izs:ize) * SQRT6*0.75_dp*k_T(ia)
+! +tau_q(ia) * sinz(izs:ize) * SQRT6*0.75_xp*k_T(ia)
! END SELECT
! CASE(1) ! j = 1
! SELECT CASE (ip-1)
! CASE(0) ! Na01 term
! moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel) = moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel)&
-! -tau_q(ia) * sinz(izs:ize) * (k_N(ia) + 3.5_dp*k_T(ia))
+! -tau_q(ia) * sinz(izs:ize) * (k_N(ia) + 3.5_xp*k_T(ia))
! CASE(2) ! Na21 term
! moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel) = moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel)&
! -tau_q(ia) * sinz(izs:ize) * SQRT2*k_T(ia)
@@ -222,7 +222,7 @@ END SUBROUTINE compute_moments_eq_rhs
! SELECT CASE (ip-1)
! CASE(0) ! Na02 term
! moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel) = moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel)&
-! +tau_q(ia) * sinz(izs:ize) * 2._dp*k_T(ia)
+! +tau_q(ia) * sinz(izs:ize) * 2._xp*k_T(ia)
! END SELECT
! END SELECT
! ENDDO
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index 37b3347a..73352402 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -10,7 +10,7 @@ MODULE nonlinear
local_nkx_ptr,kxarray, AA_x, inv_Nx,&
local_nz,ngz,zarray,nzgrid
USE model, ONLY : LINEARITY, CLOS, NL_CLOS, EM
- USE prec_const, ONLY : dp
+ USE prec_const, ONLY : xp
USE species, ONLY : sqrt_tau_o_sigma
USE time_integration, ONLY : updatetlevel
use, intrinsic :: iso_c_binding
@@ -19,12 +19,12 @@ MODULE nonlinear
INCLUDE 'fftw3-mpi.f03'
- COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: F_cmpx, G_cmpx
- COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: Fx_cmpx, Gy_cmpx
- COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: Fy_cmpx, Gx_cmpx, F_conv_G
+ COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: F_cmpx, G_cmpx
+ COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: Fx_cmpx, Gy_cmpx
+ COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: Fy_cmpx, Gx_cmpx, F_conv_G
INTEGER :: in, is, p_int, j_int, n_int
INTEGER :: nmax, smax
- REAL(dp):: sqrt_p, sqrt_pp1
+ REAL(xp):: sqrt_p, sqrt_pp1
PUBLIC :: compute_Sapj, nonlinear_init
CONTAINS
@@ -52,7 +52,7 @@ SUBROUTINE compute_Sapj
CASE ('nonlinear')
CALL compute_nonlinear
CASE ('linear')
- Sapj = 0._dp
+ Sapj = 0._xp
CASE DEFAULT
ERROR STOP '>> ERROR << Linearity not recognized '
END SELECT
@@ -69,8 +69,8 @@ SUBROUTINE compute_nonlinear
DO ip = 1,local_np ! Loop over Hermite moments
ipi = ip + ngp/2
p_int = parray(ipi)
- sqrt_p = SQRT(REAL(p_int,dp))
- sqrt_pp1 = SQRT(REAL(p_int,dp) + 1._dp)
+ sqrt_p = SQRT(REAL(p_int,xp))
+ sqrt_pp1 = SQRT(REAL(p_int,xp) + 1._xp)
eo = min(nzgrid,MODULO(parray(ip),2)+1)
DO ia = 1,local_na
IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmax)) THEN !compute for every moments except for closure 1
@@ -82,14 +82,14 @@ SUBROUTINE compute_nonlinear
ELSE
nmax = min(NL_CLOS,Jmax-j_int)
ENDIF
- bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
+ bracket_sum_r = 0._xp ! initialize sum over real nonlinear term
DO in = 1,nmax+1 ! Loop over laguerre for the sum
ini = in+ngj/2
!-----------!! ELECTROSTATIC CONTRIBUTION
! First convolution terms
F_cmpx(:,:) = phi(:,:,izi) * kernel(ia,ini,:,:,izi,eo)
! Second convolution terms
- G_cmpx = 0._dp ! initialization of the sum
+ G_cmpx = 0._xp ! initialization of the sum
smax = MIN( (in-1)+(ij-1), Jmax );
DO is = 1, smax+1 ! sum truncation on number of moments
isi = is + ngj/2
@@ -103,7 +103,7 @@ SUBROUTINE compute_nonlinear
! First convolution terms
F_cmpx(:,:) = -sqrt_tau_o_sigma(ia) * psi(:,:,izi) * kernel(ia,ini,:,:,izi,eo)
! Second convolution terms
- G_cmpx = 0._dp ! initialization of the sum
+ G_cmpx = 0._xp ! initialization of the sum
smax = MIN( (in-1)+(ij-1), Jmax );
DO is = 1, smax+1 ! sum truncation on number of moments
isi = is + ngj/2
@@ -124,7 +124,7 @@ SUBROUTINE compute_nonlinear
ENDDO
ENDDO
ELSE
- Sapj(ia,ip,ij,:,:,iz) = 0._dp
+ Sapj(ia,ip,ij,:,:,iz) = 0._xp
ENDIF
ENDDO
ENDDO
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 9b597509..9111ea3f 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -49,20 +49,20 @@ END SUBROUTINE build_dnjs_table
SUBROUTINE build_dv4Hp_table
USE array, ONLY: dv4_Hp_coeff
USE grid, ONLY: pmax
- USE prec_const, ONLY: dp, PI
+ USE prec_const, ONLY: xp, PI
IMPLICIT NONE
INTEGER :: p_
DO p_ = -2,pmax
if (p_ < 4) THEN
- dv4_Hp_coeff(p_) = 0._dp
+ dv4_Hp_coeff(p_) = 0._xp
ELSE
- dv4_Hp_coeff(p_) = 4_dp*SQRT(REAL((p_-3)*(p_-2)*(p_-1)*p_,dp))
+ dv4_Hp_coeff(p_) = 4_xp*SQRT(REAL((p_-3)*(p_-2)*(p_-1)*p_,xp))
ENDIF
ENDDO
!we scale it w.r.t. to the max degree since
!D_4^{v}\sim (\Delta v/2)^4 and \Delta v \sim 2pi/kvpar = pi/\sqrt{2P}
- ! dv4_Hp_coeff = dv4_Hp_coeff*(1._dp/2._dp/SQRT(REAL(pmax,dp)))**4
- dv4_Hp_coeff = dv4_Hp_coeff*(PI/2._dp/SQRT(2._dp*REAL(pmax,dp)))**4
+ ! dv4_Hp_coeff = dv4_Hp_coeff*(1._xp/2._xp/SQRT(REAL(pmax,xp)))**4
+ dv4_Hp_coeff = dv4_Hp_coeff*(PI/2._xp/SQRT(2._xp*REAL(pmax,xp)))**4
END SUBROUTINE build_dv4Hp_table
!******************************************************************************!
@@ -75,10 +75,10 @@ SUBROUTINE evaluate_kernels
USE grid, ONLY : local_na, local_nj,ngj, local_nkx, local_nky, local_nz, ngz, jarray, kparray,&
nzgrid
USE species, ONLY : sigma2_tau_o2
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
IMPLICIT NONE
INTEGER :: j_int, ia, eo, ikx, iky, iz, ij
- REAL(dp) :: j_dp, y_, factj
+ REAL(xp) :: j_xp, y_, factj
DO ia = 1,local_na
DO eo = 1,nzgrid
@@ -87,12 +87,12 @@ DO ia = 1,local_na
DO iz = 1,local_nz + ngz
DO ij = 1,local_nj + ngj
j_int = jarray(ij)
- j_dp = REAL(j_int,dp)
+ j_xp = REAL(j_int,xp)
y_ = sigma2_tau_o2(ia) * kparray(iky,ikx,iz,eo)**2
IF(j_int .LT. 0) THEN !ghosts values
- kernel(ia,ij,iky,ikx,iz,eo) = 0._dp
+ kernel(ia,ij,iky,ikx,iz,eo) = 0._xp
ELSE
- factj = REAL(GAMMA(j_dp+1._dp),dp)
+ factj = REAL(GAMMA(j_xp+1._xp),xp)
kernel(ia,ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/factj
ENDIF
ENDDO
@@ -102,17 +102,17 @@ DO ia = 1,local_na
ENDDO
! !! Correction term for the evaluation of the heat flux
! HF_phi_correction_operator(:,:,:) = &
- ! 2._dp * Kernel(ia,1,:,:,:,1) &
- ! -1._dp * Kernel(ia,2,:,:,:,1)
+ ! 2._xp * Kernel(ia,1,:,:,:,1) &
+ ! -1._xp * Kernel(ia,2,:,:,:,1)
!
! DO ij = 1,local_Nj
! j_int = jarray(ij)
- ! j_dp = REAL(j_int,dp)
+ ! j_xp = REAL(j_int,xp)
! HF_phi_correction_operator(:,:,:) = HF_phi_correction_operator(:,:,:) &
! - Kernel(ia,ij,:,:,:,1) * (&
- ! 2._dp*(j_dp+1.5_dp) * Kernel(ia,ij ,:,:,:,1) &
- ! - (j_dp+1.0_dp) * Kernel(ia,ij+1,:,:,:,1) &
- ! - j_dp * Kernel(ia,ij-1,:,:,:,1))
+ ! 2._xp*(j_xp+1.5_xp) * Kernel(ia,ij ,:,:,:,1) &
+ ! - (j_xp+1.0_xp) * Kernel(ia,ij+1,:,:,:,1) &
+ ! - j_xp * Kernel(ia,ij-1,:,:,:,1))
! ENDDO
ENDDO
END SUBROUTINE evaluate_kernels
@@ -135,26 +135,26 @@ SUBROUTINE evaluate_poisson_op
kxarray, kyarray, local_nj, ngj, ngz, ieven
USE species, ONLY : q2_tau
USE model, ONLY : ADIAB_E
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
IMPLICIT NONE
- REAL(dp) :: pol_tot, operator, operator_ion ! (Z^2/tau (1-sum_n kernel_na^2))
+ REAL(xp) :: pol_tot, operator, operator_ion ! (Z^2/tau (1-sum_n kernel_na^2))
INTEGER :: in,ikx,iky,iz,ia
- REAL(dp) :: sumker ! (Z_a^2/tau_a (1-sum_n kernel_na^2))
+ REAL(xp) :: sumker ! (Z_a^2/tau_a (1-sum_n kernel_na^2))
! This term has no staggered grid dependence. It is evalued for the
! even z grid since poisson uses p=0 moments and phi only.
kxloop: DO ikx = 1,local_nkx
kyloop: DO iky = 1,local_nky
zloop: DO iz = 1,local_nz
- IF( (kxarray(ikx).EQ.0._dp) .AND. (kyarray(iky).EQ.0._dp) ) THEN
- inv_poisson_op(iky, ikx, iz) = 0._dp
- inv_pol_ion (iky, ikx, iz) = 0._dp
+ IF( (kxarray(ikx).EQ.0._xp) .AND. (kyarray(iky).EQ.0._xp) ) THEN
+ inv_poisson_op(iky, ikx, iz) = 0._xp
+ inv_pol_ion (iky, ikx, iz) = 0._xp
ELSE
! loop over n only up to the max polynomial degree
- pol_tot = 0._dp ! total polarisation term
+ pol_tot = 0._xp ! total polarisation term
a:DO ia = 1,local_na ! sum over species
! ia = 1
- sumker = 0._dp ! sum of ion polarisation term
+ sumker = 0._xp ! sum of ion polarisation term
DO in=1,local_nj
sumker = sumker + q2_tau(ia)*kernel(ia,in+ngj/2,iky,ikx,iz+ngz/2,ieven)**2 ! ... sum recursively ...
END DO
@@ -162,11 +162,11 @@ ELSE
ENDDO a
operator_ion = pol_tot
IF(ADIAB_E) THEN ! Adiabatic electron model
- pol_tot = pol_tot + 1._dp
+ pol_tot = pol_tot + 1._xp
ENDIF
operator = pol_tot
- inv_poisson_op(iky, ikx, iz) = 1._dp/pol_tot
- inv_pol_ion (iky, ikx, iz) = 1._dp/operator_ion
+ inv_poisson_op(iky, ikx, iz) = 1._xp/pol_tot
+ inv_pol_ion (iky, ikx, iz) = 1._xp/operator_ion
ENDIF
END DO zloop
END DO kyloop
@@ -178,16 +178,16 @@ END SUBROUTINE evaluate_poisson_op
!!!!!!! Evaluate inverse polarisation operator for Poisson equation
!******************************************************************************!
SUBROUTINE evaluate_ampere_op
- USE prec_const, ONLY : dp
+ USE prec_const, ONLY : xp
USE array, ONLY : kernel, inv_ampere_op
USE grid, ONLY : local_na, local_nkx, local_nky, local_nz, ngz, total_nj, ngj,&
kparray, kxarray, kyarray, SOLVE_AMPERE, iodd
USE model, ONLY : beta
USE species, ONLY : q, sigma
USE geometry, ONLY : hatB
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
IMPLICIT NONE
- REAL(dp) :: sum_jpol, kperp2, operator ! (Z^2/tau (1-sum_n kernel_na^2))
+ REAL(xp) :: sum_jpol, kperp2, operator ! (Z^2/tau (1-sum_n kernel_na^2))
INTEGER :: in,ikx,iky,iz,ia
! We do not solve Ampere if beta = 0 to spare waste of ressources
IF(SOLVE_AMPERE) THEN
@@ -195,18 +195,18 @@ SUBROUTINE evaluate_ampere_op
y:DO iky = 1,local_nky
z:DO iz = 1,local_nz
kperp2 = kparray(iky,ikx,iz+ngz/2,iodd)**2
- IF( (kxarray(ikx).EQ.0._dp) .AND. (kyarray(iky).EQ.0._dp) ) THEN
- inv_ampere_op(iky, ikx, iz) = 0._dp
+ IF( (kxarray(ikx).EQ.0._xp) .AND. (kyarray(iky).EQ.0._xp) ) THEN
+ inv_ampere_op(iky, ikx, iz) = 0._xp
ELSE
- sum_jpol = 0._dp
+ sum_jpol = 0._xp
a:DO ia = 1,local_na
! loop over n only up to the max polynomial degree
j:DO in=1,total_nj
sum_jpol = sum_jpol + q(ia)**2/(sigma(ia)**2)*kernel(ia,in+ngj/2,iky,ikx,iz+ngz/2,iodd)**2 ! ... sum recursively ...
END DO j
END DO a
- operator = 2._dp*kperp2*hatB(iz+ngz/2,iodd)**2 + beta*sum_jpol
- inv_ampere_op(iky, ikx, iz) = 1._dp/operator
+ operator = 2._xp*kperp2*hatB(iz+ngz/2,iodd)**2 + beta*sum_jpol
+ inv_ampere_op(iky, ikx, iz) = 1._xp/operator
ENDIF
END DO z
END DO y
@@ -226,49 +226,49 @@ SUBROUTINE compute_lin_coeff
xpsij, xpsijp1, xpsijm1
USE species, ONLY: k_T, k_N, tau, q, sqrt_tau_o_sigma
USE model, ONLY: k_cB, k_gB
- USE prec_const, ONLY: dp, SQRT2, SQRT3
+ USE prec_const, ONLY: xp, SQRT2, SQRT3
USE grid, ONLY: parray, jarray, local_na, local_np, local_nj, ngj, ngp
INTEGER :: ia,ip,ij,p_int, j_int ! polynom. dagrees
- REAL(dp) :: p_dp, j_dp
+ REAL(xp) :: p_xp, j_xp
!! linear coefficients for moment RHS !!!!!!!!!!
DO ia = 1,local_na
DO ip = 1,local_np
p_int= parray(ip+ngp/2) ! Hermite degree
- p_dp = REAL(p_int,dp) ! REAL of Hermite degree
+ p_xp = REAL(p_int,xp) ! REAL of Hermite degree
DO ij = 1,local_nj
j_int= jarray(ij+ngj/2) ! Laguerre degree
- j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
+ j_xp = REAL(j_int,xp) ! REAL of Laguerre degree
! All Napj terms
- xnapj(ia,ip,ij) = tau(ia)/q(ia)*(k_cB*(2._dp*p_dp + 1._dp) &
- +k_gB*(2._dp*j_dp + 1._dp))
+ xnapj(ia,ip,ij) = tau(ia)/q(ia)*(k_cB*(2._xp*p_xp + 1._xp) &
+ +k_gB*(2._xp*j_xp + 1._xp))
! Mirror force terms
- ynapp1j (ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_dp+1._dp)*SQRT(p_dp+1._dp)
- ynapm1j (ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_dp+1._dp)*SQRT(p_dp)
- ynapp1jm1(ia,ip,ij) = +sqrt_tau_o_sigma(ia) * j_dp*SQRT(p_dp+1._dp)
- ynapm1jm1(ia,ip,ij) = +sqrt_tau_o_sigma(ia) * j_dp*SQRT(p_dp)
+ ynapp1j (ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_xp+1._xp)*SQRT(p_xp+1._xp)
+ ynapm1j (ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_xp+1._xp)*SQRT(p_xp)
+ ynapp1jm1(ia,ip,ij) = +sqrt_tau_o_sigma(ia) * j_xp*SQRT(p_xp+1._xp)
+ ynapm1jm1(ia,ip,ij) = +sqrt_tau_o_sigma(ia) * j_xp*SQRT(p_xp)
! Trapping terms
- zNapm1j (ia,ip,ij) = +sqrt_tau_o_sigma(ia) *(2._dp*j_dp+1._dp)*SQRT(p_dp)
- zNapm1jp1(ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_dp+1._dp)*SQRT(p_dp)
- zNapm1jm1(ia,ip,ij) = -sqrt_tau_o_sigma(ia) * j_dp*SQRT(p_dp)
+ zNapm1j (ia,ip,ij) = +sqrt_tau_o_sigma(ia) *(2._xp*j_xp+1._xp)*SQRT(p_xp)
+ zNapm1jp1(ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_xp+1._xp)*SQRT(p_xp)
+ zNapm1jm1(ia,ip,ij) = -sqrt_tau_o_sigma(ia) * j_xp*SQRT(p_xp)
ENDDO
ENDDO
DO ip = 1,local_np
p_int= parray(ip+ngp/2) ! Hermite degree
- p_dp = REAL(p_int,dp) ! REAL of Hermite degree
+ p_xp = REAL(p_int,xp) ! REAL of Hermite degree
! Landau damping coefficients (ddz napj term)
- xnapp1j(ia,ip) = sqrt_tau_o_sigma(ia) * SQRT(p_dp+1._dp)
- xnapm1j(ia,ip) = sqrt_tau_o_sigma(ia) * SQRT(p_dp)
+ xnapp1j(ia,ip) = sqrt_tau_o_sigma(ia) * SQRT(p_xp+1._xp)
+ xnapm1j(ia,ip) = sqrt_tau_o_sigma(ia) * SQRT(p_xp)
! Magnetic curvature term
- xnapp2j(ia,ip) = tau(ia)/q(ia) * k_cB * SQRT((p_dp+1._dp)*(p_dp + 2._dp))
- xnapm2j(ia,ip) = tau(ia)/q(ia) * k_cB * SQRT( p_dp *(p_dp - 1._dp))
+ xnapp2j(ia,ip) = tau(ia)/q(ia) * k_cB * SQRT((p_xp+1._xp)*(p_xp + 2._xp))
+ xnapm2j(ia,ip) = tau(ia)/q(ia) * k_cB * SQRT( p_xp *(p_xp - 1._xp))
ENDDO
DO ij = 1,local_nj
j_int= jarray(ij+ngj/2) ! Laguerre degree
- j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
+ j_xp = REAL(j_int,xp) ! REAL of Laguerre degree
! Magnetic gradient term
- xnapjp1(ia,ij) = -tau(ia)/q(ia) * k_gB * (j_dp + 1._dp)
- xnapjm1(ia,ij) = -tau(ia)/q(ia) * k_gB * j_dp
+ xnapjp1(ia,ij) = -tau(ia)/q(ia) * k_gB * (j_xp + 1._xp)
+ xnapjm1(ia,ij) = -tau(ia)/q(ia) * k_gB * j_xp
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! ES linear coefficients for moment RHS !!!!!!!!!!
@@ -276,18 +276,18 @@ SUBROUTINE compute_lin_coeff
p_int= parray(ip+ngp/2) ! Hermite degree
DO ij = 1,local_nj
j_int= jarray(ij+ngj/2) ! REALof Laguerre degree
- j_dp = REAL(j_int,dp) ! REALof Laguerre degree
+ j_xp = REAL(j_int,xp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 0) THEN ! kronecker p0
- xphij (ia,ip,ij) = +k_N(ia) + 2._dp*j_dp*k_T(ia)
- xphijp1(ia,ip,ij) = -k_T(ia)*(j_dp+1._dp)
- xphijm1(ia,ip,ij) = -k_T(ia)* j_dp
+ xphij (ia,ip,ij) = +k_N(ia) + 2._xp*j_xp*k_T(ia)
+ xphijp1(ia,ip,ij) = -k_T(ia)*(j_xp+1._xp)
+ xphijm1(ia,ip,ij) = -k_T(ia)* j_xp
ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
xphij(ia,ip,ij) = +k_T(ia)/SQRT2
- xphijp1(ia,ip,ij) = 0._dp; xphijm1(ia,ip,ij) = 0._dp;
+ xphijp1(ia,ip,ij) = 0._xp; xphijm1(ia,ip,ij) = 0._xp;
ELSE
- xphij (ia,ip,ij) = 0._dp; xphijp1(ia,ip,ij) = 0._dp
- xphijm1(ia,ip,ij) = 0._dp;
+ xphij (ia,ip,ij) = 0._xp; xphijp1(ia,ip,ij) = 0._xp
+ xphijm1(ia,ip,ij) = 0._xp;
ENDIF
ENDDO
ENDDO
@@ -297,17 +297,17 @@ SUBROUTINE compute_lin_coeff
p_int= parray(ip+ngp/2) ! Hermite degree
DO ij = 1,local_nj
j_int= jarray(ij+ngj/2) ! REALof Laguerre degree
- j_dp = REAL(j_int,dp) ! REALof Laguerre degree
+ j_xp = REAL(j_int,xp) ! REALof Laguerre degree
IF (p_int .EQ. 1) THEN ! kronecker p1
- xpsij (ia,ip,ij) = +(k_N(ia) + (2._dp*j_dp+1._dp)*k_T(ia))* sqrt_tau_o_sigma(ia)
- xpsijp1(ia,ip,ij) = - k_T(ia)*(j_dp+1._dp) * sqrt_tau_o_sigma(ia)
- xpsijm1(ia,ip,ij) = - k_T(ia)* j_dp * sqrt_tau_o_sigma(ia)
+ xpsij (ia,ip,ij) = +(k_N(ia) + (2._xp*j_xp+1._xp)*k_T(ia))* sqrt_tau_o_sigma(ia)
+ xpsijp1(ia,ip,ij) = - k_T(ia)*(j_xp+1._xp) * sqrt_tau_o_sigma(ia)
+ xpsijm1(ia,ip,ij) = - k_T(ia)* j_xp * sqrt_tau_o_sigma(ia)
ELSE IF (p_int .EQ. 3) THEN ! kronecker p3
xpsij (ia,ip,ij) = + k_T(ia)*SQRT3/SQRT2 * sqrt_tau_o_sigma(ia)
- xpsijp1(ia,ip,ij) = 0._dp; xpsijm1(ia,ip,ij) = 0._dp;
+ xpsijp1(ia,ip,ij) = 0._xp; xpsijm1(ia,ip,ij) = 0._xp;
ELSE
- xpsij (ia,ip,ij) = 0._dp; xpsijp1(ia,ip,ij) = 0._dp
- xpsijm1(ia,ip,ij) = 0._dp;
+ xpsij (ia,ip,ij) = 0._xp; xpsijp1(ia,ip,ij) = 0._xp
+ xpsijm1(ia,ip,ij) = 0._xp;
ENDIF
ENDDO
ENDDO
diff --git a/src/parallel_mod.F90 b/src/parallel_mod.F90
index e4469d0f..0dbb02cb 100644
--- a/src/parallel_mod.F90
+++ b/src/parallel_mod.F90
@@ -1,5 +1,5 @@
MODULE parallel
- use prec_const, ONLY : dp
+ use prec_const, ONLY : xp
USE mpi
IMPLICIT NONE
! Auxiliary variables
@@ -239,12 +239,12 @@ CONTAINS
SUBROUTINE gather_xyz(field_loc,field_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot)
IMPLICIT NONE
INTEGER, INTENT(IN) :: nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot
- COMPLEX(dp), DIMENSION(nky_loc,nkx_tot,nz_loc), INTENT(IN) :: field_loc
- COMPLEX(dp), DIMENSION(nky_tot,nkx_tot,nz_tot), INTENT(OUT) :: field_tot
- COMPLEX(dp), DIMENSION(nky_tot,nz_loc) :: buffer_yt_zl !full y, local z
- COMPLEX(dp), DIMENSION(nky_tot,nz_tot) :: buffer_yt_zt !full y, full z
- COMPLEX(dp), DIMENSION(nky_loc):: buffer_yl_zc !local y, constant z
- COMPLEX(dp), DIMENSION(nky_tot):: buffer_yt_zc !full y, constant z
+ COMPLEX(xp), DIMENSION(nky_loc,nkx_tot,nz_loc), INTENT(IN) :: field_loc
+ COMPLEX(xp), DIMENSION(nky_tot,nkx_tot,nz_tot), INTENT(OUT) :: field_tot
+ COMPLEX(xp), DIMENSION(nky_tot,nz_loc) :: buffer_yt_zl !full y, local z
+ COMPLEX(xp), DIMENSION(nky_tot,nz_tot) :: buffer_yt_zt !full y, full z
+ COMPLEX(xp), DIMENSION(nky_loc):: buffer_yl_zc !local y, constant z
+ COMPLEX(xp), DIMENSION(nky_tot):: buffer_yt_zc !full y, constant z
INTEGER :: snd_y, snd_z, root_p, root_z, root_ky, ix, iz
snd_y = nky_loc ! Number of points to send along y (per z)
@@ -277,12 +277,12 @@ CONTAINS
SUBROUTINE gather_pjz(field_loc,field_tot,np_loc,np_tot,nj_tot,nz_loc,nz_tot)
IMPLICIT NONE
INTEGER, INTENT(IN) :: np_loc,np_tot, nj_tot, nz_loc,nz_tot
- COMPLEX(dp), DIMENSION(np_loc,nj_tot,nz_loc), INTENT(IN) :: field_loc
- COMPLEX(dp), DIMENSION(np_tot,nj_tot,nz_tot), INTENT(OUT) :: field_tot
- COMPLEX(dp), DIMENSION(np_tot,nz_loc) :: buffer_pt_zl !full p, local z
- COMPLEX(dp), DIMENSION(np_tot,nz_tot) :: buffer_pt_zt !full p, full z
- COMPLEX(dp), DIMENSION(np_loc) :: buffer_pl_zc !local p, constant z
- COMPLEX(dp), DIMENSION(np_tot) :: buffer_pt_zc !full p, constant z
+ COMPLEX(xp), DIMENSION(np_loc,nj_tot,nz_loc), INTENT(IN) :: field_loc
+ COMPLEX(xp), DIMENSION(np_tot,nj_tot,nz_tot), INTENT(OUT) :: field_tot
+ COMPLEX(xp), DIMENSION(np_tot,nz_loc) :: buffer_pt_zl !full p, local z
+ COMPLEX(xp), DIMENSION(np_tot,nz_tot) :: buffer_pt_zt !full p, full z
+ COMPLEX(xp), DIMENSION(np_loc) :: buffer_pl_zc !local p, constant z
+ COMPLEX(xp), DIMENSION(np_tot) :: buffer_pt_zc !full p, constant z
INTEGER :: snd_p, snd_z, root_p, root_z, root_ky, ij, iz
snd_p = np_loc ! Number of points to send along p (per z)
@@ -317,14 +317,14 @@ CONTAINS
SUBROUTINE gather_pjxyz(field_loc,field_tot,na_tot,np_loc,np_tot,nj_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot)
IMPLICIT NONE
INTEGER, INTENT(IN) :: na_tot,np_loc,np_tot,nj_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot
- COMPLEX(dp), DIMENSION(np_loc,nj_tot,nky_loc,nkx_tot,nz_loc), INTENT(IN) :: field_loc
- COMPLEX(dp), DIMENSION(na_tot,np_tot,nj_tot,nky_tot,nkx_tot,nz_tot), INTENT(OUT) :: field_tot
- COMPLEX(dp), DIMENSION(np_tot,nky_tot,nz_loc) :: buffer_pt_yt_zl ! full p, full y, local z
- COMPLEX(dp), DIMENSION(np_tot,nky_tot,nz_tot) :: buffer_pt_yt_zt ! full p, full y, full z
- COMPLEX(dp), DIMENSION(np_tot,nky_loc) :: buffer_pt_yl_zc ! full p, local y, constant z
- COMPLEX(dp), DIMENSION(np_tot,nky_tot) :: buffer_pt_yt_zc ! full p, full y, constant z
- COMPLEX(dp), DIMENSION(np_loc) :: buffer_pl_cy_zc !local p, constant y, constant z
- COMPLEX(dp), DIMENSION(np_tot) :: buffer_pt_cy_zc ! full p, constant y, constant z
+ COMPLEX(xp), DIMENSION(np_loc,nj_tot,nky_loc,nkx_tot,nz_loc), INTENT(IN) :: field_loc
+ COMPLEX(xp), DIMENSION(na_tot,np_tot,nj_tot,nky_tot,nkx_tot,nz_tot), INTENT(OUT) :: field_tot
+ COMPLEX(xp), DIMENSION(np_tot,nky_tot,nz_loc) :: buffer_pt_yt_zl ! full p, full y, local z
+ COMPLEX(xp), DIMENSION(np_tot,nky_tot,nz_tot) :: buffer_pt_yt_zt ! full p, full y, full z
+ COMPLEX(xp), DIMENSION(np_tot,nky_loc) :: buffer_pt_yl_zc ! full p, local y, constant z
+ COMPLEX(xp), DIMENSION(np_tot,nky_tot) :: buffer_pt_yt_zc ! full p, full y, constant z
+ COMPLEX(xp), DIMENSION(np_loc) :: buffer_pl_cy_zc !local p, constant y, constant z
+ COMPLEX(xp), DIMENSION(np_tot) :: buffer_pt_cy_zc ! full p, constant y, constant z
INTEGER :: snd_p, snd_y, snd_z, root_p, root_z, root_ky, iy, ix, iz, ij, ia
snd_p = np_loc ! Number of points to send along p (per z,y)
snd_y = np_tot*nky_loc ! Number of points to send along y (per z, full p)
@@ -368,8 +368,8 @@ CONTAINS
SUBROUTINE manual_3D_bcast(field_,n1,n2,n3)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n1,n2,n3
- COMPLEX(dp), DIMENSION(n1,n2,n3), INTENT(INOUT) :: field_
- COMPLEX(dp) :: buffer(n1,n2,n3)
+ COMPLEX(xp), DIMENSION(n1,n2,n3), INTENT(INOUT) :: field_
+ COMPLEX(xp) :: buffer(n1,n2,n3)
INTEGER :: i_, root, world_rank, world_size, count, i1,i2,i3
root = 0;
count = n1*n2*n3;
@@ -410,8 +410,8 @@ CONTAINS
!!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!!
SUBROUTINE manual_0D_bcast(v)
IMPLICIT NONE
- COMPLEX(dp), INTENT(INOUT) :: v
- COMPLEX(dp) :: buffer
+ COMPLEX(xp), INTENT(INOUT) :: v
+ COMPLEX(xp) :: buffer
INTEGER :: i_, root, world_rank, world_size, count
root = 0;
count = 1;
@@ -441,7 +441,7 @@ CONTAINS
SUBROUTINE exchange_ghosts_1D(f,nbr_L,nbr_R,np,ng)
IMPLICIT NONE
INTEGER, INTENT(IN) :: np,ng, nbr_L, nbr_R
- COMPLEX(dp), DIMENSION(np+ng), INTENT(INOUT) :: f
+ COMPLEX(xp), DIMENSION(np+ng), INTENT(INOUT) :: f
INTEGER :: ierr, first, last, ig
first = 1 + ng/2
last = np + ng/2
diff --git a/src/prec_const_mod.F90 b/src/prec_const_mod.F90
index ba8a8dd7..af0ac57c 100644
--- a/src/prec_const_mod.F90
+++ b/src/prec_const_mod.F90
@@ -6,43 +6,50 @@ MODULE prec_const
stdin=>input_unit, &
stdout=>output_unit, &
stderr=>error_unit
- ! Precision for real and complex
- INTEGER, PARAMETER :: sp = REAL32 !Single precision, should not be used
- INTEGER, PARAMETER :: dp = REAL64 !real(dp), enforced through the code
-
- INTEGER, private :: dp_r, dp_p !Range and Aprecision of doubles
- INTEGER, private :: sp_r, sp_p !Range and precision of singles
-
-
- INTEGER, private :: MPI_SP !Single precision for MPI
- INTEGER, private :: MPI_DP !Double precision in MPI
- INTEGER, private :: MPI_SUM_DP !Sum reduction operation for DP datatype
- INTEGER, private :: MPI_MAX_DP !Max reduction operation for DP datatype
- INTEGER, private :: MPI_MIN_DP !Min reduction operation for DP datatype
+ use, intrinsic :: iso_c_binding
+
+ ! Define single and double precision
+ ! INTEGER, PARAMETER :: sp = REAL32 !Single precision
+ ! INTEGER, PARAMETER :: dp = REAL64 !Double precision
+ ! INTEGER, private :: dp_r, dp_p !Range and Aprecision of doubles
+ ! INTEGER, private :: sp_r, sp_p !Range and precision of singles
+ ! INTEGER, private :: MPI_SP !Single precision for MPI
+ ! INTEGER, private :: MPI_DP !Double precision in MPI
+ ! INTEGER, private :: MPI_SUM_DP !Sum reduction operation for DP datatype
+ ! INTEGER, private :: MPI_MAX_DP !Max reduction operation for DP datatype
+ ! INTEGER, private :: MPI_MIN_DP !Min reduction operation for DP datatype
+
+ ! Define a generic precision parameter for the entire program
+ INTEGER, PARAMETER :: xp = REAL64!real(xp), enforced through the code
+ INTEGER, private :: xp_r, xp_p !Range and precision of single
+ INTEGER, private :: MPI_XP !Double precision in MPI
+ INTEGER, private :: MPI_SUM_XP !Sum reduction operation for xp datatype
+ INTEGER, private :: MPI_MAX_XP !Max reduction operation for xp datatype
+ INTEGER, private :: MPI_MIN_XP !Min reduction operation for xp datatype
! Some useful constants, to avoid recomputing them too often
- REAL(dp), PARAMETER :: PI=3.141592653589793238462643383279502884197_dp
- REAL(dp), PARAMETER :: PIO2=1.57079632679489661923132169163975144209858_dp
- REAL(dp), PARAMETER :: TWOPI=6.283185307179586476925286766559005768394_dp
- REAL(dp), PARAMETER :: ONEOPI=0.3183098861837906912164442019275156781077_dp
- REAL(dp), PARAMETER :: SQRT2=1.41421356237309504880168872420969807856967_dp
- REAL(dp), PARAMETER :: SQRT6=SQRT(6._dp)
- REAL(dp), PARAMETER :: INVSQRT2=0.7071067811865475244008443621048490392848359377_dp
- REAL(dp), PARAMETER :: SQRT3=1.73205080756887729352744634150587236694281_dp
- REAL(dp), PARAMETER :: onetwelfth=0.08333333333333333333333333333333333333333333333_dp
- REAL(dp), PARAMETER :: onetwentyfourth=0.04166666666666666666666666666666666666666666667_dp
- REAL(dp), PARAMETER :: onethird=0.33333333333333333333333333333333333333333333333_dp
- REAL(dp), PARAMETER :: twothird=0.66666666666666666666666666666666666666666666666_dp
- REAL(dp), PARAMETER :: onesixth=0.1666666666666666666666666666666666666666666667_dp
- REAL(dp), PARAMETER :: fivesixths=0.8333333333333333333333333333333333333333333333_dp
- REAL(dp), PARAMETER :: sevensixths=1.1666666666666666666666666666666666666666666667_dp
- REAL(dp), PARAMETER :: elevensixths=1.833333333333333333333333333333333333333333333_dp
- REAL(dp), PARAMETER :: nineeighths=1.125_dp
- REAL(dp), PARAMETER :: onesixteenth=0.0625_dp
- REAL(dp), PARAMETER :: ninesixteenths=0.5625_dp
- REAL(dp), PARAMETER :: thirteentwelfths = 1.083333333333333333333333333333333333333333333_dp
- COMPLEX(dp), PARAMETER :: imagu = (0._dp,1._dp)
+ REAL(xp), PARAMETER :: PI=3.141592653589793238462643383279502884197_xp
+ REAL(xp), PARAMETER :: PIO2=1.57079632679489661923132169163975144209858_xp
+ REAL(xp), PARAMETER :: TWOPI=6.283185307179586476925286766559005768394_xp
+ REAL(xp), PARAMETER :: ONEOPI=0.3183098861837906912164442019275156781077_xp
+ REAL(xp), PARAMETER :: SQRT2=1.41421356237309504880168872420969807856967_xp
+ REAL(xp), PARAMETER :: SQRT6=SQRT(6._xp)
+ REAL(xp), PARAMETER :: INVSQRT2=0.7071067811865475244008443621048490392848359377_xp
+ REAL(xp), PARAMETER :: SQRT3=1.73205080756887729352744634150587236694281_xp
+ REAL(xp), PARAMETER :: onetwelfth=0.08333333333333333333333333333333333333333333333_xp
+ REAL(xp), PARAMETER :: onetwentyfourth=0.04166666666666666666666666666666666666666666667_xp
+ REAL(xp), PARAMETER :: onethird=0.33333333333333333333333333333333333333333333333_xp
+ REAL(xp), PARAMETER :: twothird=0.66666666666666666666666666666666666666666666666_xp
+ REAL(xp), PARAMETER :: onesixth=0.1666666666666666666666666666666666666666666667_xp
+ REAL(xp), PARAMETER :: fivesixths=0.8333333333333333333333333333333333333333333333_xp
+ REAL(xp), PARAMETER :: sevensixths=1.1666666666666666666666666666666666666666666667_xp
+ REAL(xp), PARAMETER :: elevensixths=1.833333333333333333333333333333333333333333333_xp
+ REAL(xp), PARAMETER :: nineeighths=1.125_xp
+ REAL(xp), PARAMETER :: onesixteenth=0.0625_xp
+ REAL(xp), PARAMETER :: ninesixteenths=0.5625_xp
+ REAL(xp), PARAMETER :: thirteentwelfths = 1.083333333333333333333333333333333333333333333_xp
+ COMPLEX(xp), PARAMETER :: imagu = (0._xp,1._xp)
CONTAINS
SUBROUTINE INIT_PREC_CONST
@@ -50,21 +57,24 @@ MODULE prec_const
IMPLICIT NONE
integer :: ierr,me
- REAL(sp) :: a = 1_sp
- REAL(dp) :: b = 1_dp
-
+ ! REAL(sp) :: a = 1_sp
+ ! REAL(dp) :: b = 1_dp
!Get range and precision of ISO FORTRAN sizes
- sp_r = range(a)
- sp_p = precision(a)
-
- dp_r = range(b)
- dp_p = precision(b)
+ ! sp_r = range(a)
+ ! sp_p = precision(a)
+ ! dp_r = range(b)
+ ! dp_p = precision(b)
+
+ REAL(xp) :: c = 1_xp
+ xp_r = range(c)
+ xp_p = precision(c)
CALL mpi_comm_rank(MPI_COMM_WORLD,me,ierr)
!Create MPI datatypes that support the specific size
- CALL MPI_Type_create_f90_real(sp_p,sp_r,MPI_sp,ierr)
- CALL MPI_Type_create_f90_real(dp_p,dp_r,MPI_dp,ierr)
+ ! CALL MPI_Type_create_f90_real(sp_p,sp_r,MPI_sp,ierr)
+ ! CALL MPI_Type_create_f90_real(dp_p,dp_r,MPI_xp,ierr)
+ CALL MPI_Type_create_f90_real(xp_p,xp_r,MPI_xp,ierr)
END SUBROUTINE INIT_PREC_CONST
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 46532eaa..c7e8dbee 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -1,5 +1,5 @@
MODULE processing
- USE prec_const, ONLY: dp, imagu, SQRT2, SQRT3, onetwelfth, twothird
+ USE prec_const, ONLY: xp, imagu, SQRT2, SQRT3, onetwelfth, twothird
USE grid, ONLY: &
local_na, local_np, local_nj, local_nky, local_nkx, local_nz, Ngz,Ngj,Ngp,nzgrid, &
parray,pmax,ip0, iodd, ieven,&
@@ -21,8 +21,8 @@ MODULE processing
USE mpi
implicit none
- REAL(dp), PUBLIC, ALLOCATABLE, DIMENSION(:), PROTECTED :: pflux_x, gflux_x
- REAL(dp), PUBLIC, ALLOCATABLE, DIMENSION(:), PROTECTED :: hflux_x
+ REAL(xp), PUBLIC, ALLOCATABLE, DIMENSION(:), PROTECTED :: pflux_x, gflux_x
+ REAL(xp), PUBLIC, ALLOCATABLE, DIMENSION(:), PROTECTED :: hflux_x
INTEGER :: ierr
PUBLIC :: init_process
PUBLIC :: compute_nadiab_moments, compute_gradients_z, compute_interp_z
@@ -83,7 +83,7 @@ CONTAINS
p_int = parray(ip)
DO ia = 1,local_na
IF(p_int+2*j_int .GT. dmax) &
- nadiab_moments(ia,ip,ij,iky,ikx,iz) = 0._dp
+ nadiab_moments(ia,ip,ij,iky,ikx,iz) = 0._xp
ENDDO
ENDDO
ENDDO
@@ -97,8 +97,8 @@ CONTAINS
! SUBROUTINE compute_gradients_z
! IMPLICIT NONE
! INTEGER :: eo, p_int, j_int, ia,ip,ij,iky,ikx,iz,izi
- ! COMPLEX(dp), DIMENSION(local_nz+ngz) :: f_in
- ! COMPLEX(dp), DIMENSION(local_nz) :: f_out
+ ! COMPLEX(xp), DIMENSION(local_nz+ngz) :: f_in
+ ! COMPLEX(xp), DIMENSION(local_nz) :: f_out
! ! Compute z first derivative
! DO iz=1,local_nz+ngz
! izi = iz+ngz/2
@@ -114,11 +114,11 @@ CONTAINS
! -onetwelfth*nadiab_moments(ia,ip,ij,iky,ikx,izi-2)&
! )
! ddzND_Napj(ia,ip,ij,iky,ikx,iz) = inv_deltaz**4 *(&
- ! +1._dp*moments(ia,ip,ij,iky,ikx,izi-2,updatetlevel)&
- ! -4._dp*moments(ia,ip,ij,iky,ikx,izi-1,updatetlevel)&
- ! +6._dp*moments(ia,ip,ij,iky,ikx,izi ,updatetlevel)&
- ! -4._dp*moments(ia,ip,ij,iky,ikx,izi+1,updatetlevel)&
- ! +1._dp*moments(ia,ip,ij,iky,ikx,izi-2,updatetlevel)&
+ ! +1._xp*moments(ia,ip,ij,iky,ikx,izi-2,updatetlevel)&
+ ! -4._xp*moments(ia,ip,ij,iky,ikx,izi-1,updatetlevel)&
+ ! +6._xp*moments(ia,ip,ij,iky,ikx,izi ,updatetlevel)&
+ ! -4._xp*moments(ia,ip,ij,iky,ikx,izi+1,updatetlevel)&
+ ! +1._xp*moments(ia,ip,ij,iky,ikx,izi-2,updatetlevel)&
! )
! ENDDO
! ENDDO
@@ -131,9 +131,9 @@ CONTAINS
! ! z grid gradients
SUBROUTINE compute_gradients_z
IMPLICIT NONE
- INTEGER :: eo, p_int, j_int, ia,ip,ij,iky,ikx,iz
- COMPLEX(dp), DIMENSION(local_nz+ngz) :: f_in
- COMPLEX(dp), DIMENSION(local_nz) :: f_out
+ INTEGER :: eo, p_int, ia,ip,ij,iky,ikx,iz
+ COMPLEX(xp), DIMENSION(local_nz+ngz) :: f_in
+ COMPLEX(xp), DIMENSION(local_nz) :: f_out
DO ikx = 1,local_nkx
DO iky = 1,local_nky
DO ij = 1,local_nj+ngj
@@ -170,9 +170,9 @@ CONTAINS
! z grid interpolation
SUBROUTINE compute_interp_z
IMPLICIT NONE
- INTEGER :: eo, p_int, j_int, ia,ip,ij,iky,ikx,iz
- COMPLEX(dp), DIMENSION(local_nz+ngz) :: f_in
- COMPLEX(dp), DIMENSION(local_nz) :: f_out
+ INTEGER :: eo, ia,ip,ij,iky,ikx,iz
+ COMPLEX(xp), DIMENSION(local_nz+ngz) :: f_in
+ COMPLEX(xp), DIMENSION(local_nz) :: f_out
IF(nzgrid .GT. 1) THEN
DO ikx = 1,local_nkx
DO iky = 1,local_nky
@@ -201,14 +201,14 @@ CONTAINS
! 1D diagnostic to compute the average radial particle transport <n_a v_ExB_x>_xyz
SUBROUTINE compute_radial_transport
IMPLICIT NONE
- COMPLEX(dp) :: pflux_local, gflux_local, integral
- REAL(dp) :: buffer(2)
+ COMPLEX(xp) :: pflux_local, gflux_local, integral
+ REAL(xp) :: buffer(2)
INTEGER :: i_, root, iky, ikx, ia, iz, in, izi, ini
- COMPLEX(dp), DIMENSION(local_nz) :: integrant
+ COMPLEX(xp), DIMENSION(local_nz) :: integrant
DO ia = 1,local_na
- pflux_local = 0._dp ! particle flux
- gflux_local = 0._dp ! gyrocenter flux
- integrant = 0._dp ! auxiliary variable for z integration
+ pflux_local = 0._xp ! particle flux
+ gflux_local = 0._xp ! gyrocenter flux
+ integrant = 0._xp ! auxiliary variable for z integration
!!---------- Gyro center flux (drift kinetic) ------------
! Electrostatic part
IF(CONTAINSp0) THEN
@@ -239,10 +239,10 @@ CONTAINS
! Integrate over z
call simpson_rule_z(local_nz,zweights_SR,integrant,integral)
! Get process local gyrocenter flux with a factor two to account for the negative ky modes
- gflux_local = 2._dp*integral*iInt_Jacobian
+ gflux_local = 2._xp*integral*iInt_Jacobian
!
!!---------- Particle flux (gyrokinetic) ------------
- integrant = 0._dp ! reset auxiliary variable
+ integrant = 0._xp ! reset auxiliary variable
! Electrostatic part
IF(CONTAINSp0) THEN
DO iz = 1,local_nz ! we take interior points only
@@ -278,10 +278,10 @@ CONTAINS
! Integrate over z
call simpson_rule_z(local_nz,zweights_SR,integrant,integral)
! Get process local particle flux with a factor two to account for the negative ky modes
- pflux_local = 2._dp*integral*iInt_Jacobian
+ pflux_local = 2._xp*integral*iInt_Jacobian
!!!!---------- Sum over all processes ----------
- buffer(1) = REAL(gflux_local,dp)
- buffer(2) = REAL(pflux_local,dp)
+ buffer(1) = REAL(gflux_local,xp)
+ buffer(2) = REAL(pflux_local,xp)
root = 0
!Gather manually among the rank_p=0 processes and perform the sum
gflux_x(ia) = 0
@@ -309,13 +309,13 @@ CONTAINS
! 1D diagnostic to compute the average radial particle heatflux <T_i v_ExB_x>_xyz
SUBROUTINE compute_radial_heatflux
IMPLICIT NONE
- COMPLEX(dp) :: hflux_local, integral
- REAL(dp) :: buffer(2), n_dp
+ COMPLEX(xp) :: hflux_local, integral
+ REAL(xp) :: buffer(2), n_xp
INTEGER :: i_, root, in, ia, iky, ikx, iz, izi, ini
- COMPLEX(dp), DIMENSION(local_nz) :: integrant ! charge density q_a n_a
+ COMPLEX(xp), DIMENSION(local_nz) :: integrant ! charge density q_a n_a
DO ia = 1,local_na
- hflux_local = 0._dp ! heat flux
- integrant = 0._dp ! z integration auxiliary variable
+ hflux_local = 0._xp ! heat flux
+ integrant = 0._xp ! z integration auxiliary variable
!!----------------ELECTROSTATIC CONTRIBUTION---------------------------
IF(CONTAINSp0 .AND. CONTAINSp2) THEN
! Loop to compute gamma_kx = sum_ky sum_j -i k_z Kernel_j Na00 * phi
@@ -325,14 +325,14 @@ CONTAINS
DO iky = 1,local_nky
DO in = 1, local_nj
ini = in+ngj/2 !interior index for ghosted arrays
- n_dp = jarray(ini)
+ n_xp = jarray(ini)
integrant(iz) = integrant(iz) &
-Jacobian(izi,ieven)*tau(ia)*imagu*kyarray(iky)*phi(iky,ikx,izi)&
*kernel(ia,ini,iky,ikx,izi,ieven)*CONJG(&
- 0.5_dp*SQRT2*moments(ia,ip2,ini ,iky,ikx,izi,updatetlevel)&
- +(2._dp*n_dp + 1.5_dp)*moments(ia,ip0,ini ,iky,ikx,izi,updatetlevel)&
- -(n_dp+1._dp)*moments(ia,ip0,ini+1,iky,ikx,izi,updatetlevel)&
- -n_dp*moments(ia,ip0,ini-1,iky,ikx,izi,updatetlevel))
+ 0.5_xp*SQRT2*moments(ia,ip2,ini ,iky,ikx,izi,updatetlevel)&
+ +(2._xp*n_xp + 1.5_xp)*moments(ia,ip0,ini ,iky,ikx,izi,updatetlevel)&
+ -(n_xp+1._xp)*moments(ia,ip0,ini+1,iky,ikx,izi,updatetlevel)&
+ -n_xp*moments(ia,ip0,ini-1,iky,ikx,izi,updatetlevel))
ENDDO
ENDDO
ENDDO
@@ -348,15 +348,15 @@ CONTAINS
DO iky = 1,local_nky
DO in = 1, local_nj
ini = in + ngj/2 !interior index for ghosted arrays
- n_dp = jarray(ini)
+ n_xp = jarray(ini)
integrant(iz) = integrant(iz) &
+Jacobian(izi,iodd)*tau(ia)*sqrt_tau_o_sigma(ia)*imagu*kyarray(iky)*CONJG(psi(iky,ikx,izi))&
*kernel(ia,ini,iky,ikx,izi,iodd)*(&
- 0.5_dp*SQRT2*SQRT3*moments(ia,ip3,ini ,iky,ikx,izi,updatetlevel)&
- +1.5_dp*moments(ia,ip1,ini ,iky,ikx,izi,updatetlevel)&
- +(2._dp*n_dp+1._dp)*moments(ia,ip1,ini ,iky,ikx,izi,updatetlevel)&
- -(n_dp+1._dp)*moments(ia,ip1,ini+1,iky,ikx,izi,updatetlevel)&
- -n_dp*moments(ia,ip1,ini-1,iky,ikx,izi,updatetlevel))
+ 0.5_xp*SQRT2*SQRT3*moments(ia,ip3,ini ,iky,ikx,izi,updatetlevel)&
+ +1.5_xp*moments(ia,ip1,ini ,iky,ikx,izi,updatetlevel)&
+ +(2._xp*n_xp+1._xp)*moments(ia,ip1,ini ,iky,ikx,izi,updatetlevel)&
+ -(n_xp+1._xp)*moments(ia,ip1,ini+1,iky,ikx,izi,updatetlevel)&
+ -n_xp*moments(ia,ip1,ini-1,iky,ikx,izi,updatetlevel))
ENDDO
ENDDO
ENDDO
@@ -368,8 +368,8 @@ CONTAINS
! Integrate over z
call simpson_rule_z(local_nz,zweights_SR,integrant,integral)
! Double it for taking into account the other half plane
- hflux_local = 2._dp*integral*iInt_Jacobian
- buffer(2) = REAL(hflux_local,dp)
+ hflux_local = 2._xp*integral*iInt_Jacobian
+ buffer(2) = REAL(hflux_local,xp)
root = 0
!Gather manually among the rank_p=0 processes and perform the sum
hflux_x(ia) = 0
@@ -396,13 +396,13 @@ CONTAINS
USE array, ONLY : Napjz
USE time_integration, ONLY : updatetlevel
IMPLICIT NONE
- REAL(dp), DIMENSION(local_np,local_nj,local_nz) :: local_sum,global_sum, buffer
+ REAL(xp), DIMENSION(local_np,local_nj,local_nz) :: local_sum,global_sum, buffer
INTEGER :: i_, root, count, ia, ip, ij, iky, ikx, iz
root = 0
DO ia=1,local_na
! z-moment spectrum
! build local sum
- local_sum = 0._dp
+ local_sum = 0._xp
DO iz = 1,local_nz
DO ikx = 1,local_nkx
DO iky = 1,local_nky
@@ -418,7 +418,7 @@ CONTAINS
ENDDO
! sum reduction
buffer = local_sum
- global_sum = 0._dp
+ global_sum = 0._xp
count = local_np*local_nj*local_nz
IF (num_procs_ky .GT. 1) THEN
!! Everyone sends its local_sum to root = 0
@@ -444,7 +444,7 @@ CONTAINS
! Compute the 2D particle density for electron and ions (sum over Laguerre)
SUBROUTINE compute_density
IMPLICIT NONE
- COMPLEX(dp) :: dens_
+ COMPLEX(xp) :: dens_
INTEGER :: ia, iz, iky, ikx, ij
DO ia=1,local_na
IF ( CONTAINSp0 ) THEN
@@ -452,7 +452,7 @@ CONTAINS
DO iz = 1,local_nz
DO iky = 1,local_nky
DO ikx = 1,local_nkx
- dens_ = 0._dp
+ dens_ = 0._xp
DO ij = 1, local_nj
dens_ = dens_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven) * moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
ENDDO
@@ -467,17 +467,17 @@ CONTAINS
! Compute the 2D particle fluid perp velocity for electron and ions (sum over Laguerre)
SUBROUTINE compute_uperp
IMPLICIT NONE
- COMPLEX(dp) :: uperp_
+ COMPLEX(xp) :: uperp_
INTEGER :: ia, iz, iky, ikx, ij
DO ia=1,local_na
IF ( CONTAINSp0 ) THEN
DO iz = 1,local_nz
DO iky = 1,local_nky
DO ikx = 1,local_nkx
- uperp_ = 0._dp
+ uperp_ = 0._xp
DO ij = 1, local_nj
uperp_ = uperp_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven) *&
- 0.5_dp*(moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
+ 0.5_xp*(moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
-moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
ENDDO
uper(ia,iky,ikx,iz) = uperp_
@@ -492,13 +492,13 @@ CONTAINS
SUBROUTINE compute_upar
IMPLICIT NONE
INTEGER :: ia, iz, iky, ikx, ij
- COMPLEX(dp) :: upar_
+ COMPLEX(xp) :: upar_
DO ia=1,local_na
IF ( CONTAINSp1 ) THEN
DO iz = 1,local_nz
DO iky = 1,local_nky
DO ikx = 1,local_nkx
- upar_ = 0._dp
+ upar_ = 0._xp
DO ij = 1, local_nj
upar_ = upar_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,iodd)*moments(ia,ip1,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
ENDDO
@@ -514,8 +514,8 @@ CONTAINS
SUBROUTINE compute_tperp
USE time_integration, ONLY : updatetlevel
IMPLICIT NONE
- REAL(dp) :: j_dp
- COMPLEX(dp) :: Tperp_
+ REAL(xp) :: j_xp
+ COMPLEX(xp) :: Tperp_
INTEGER :: ia, iz, iky, ikx, ij
DO ia=1,local_na
IF ( CONTAINSp0 .AND. CONTAINSp2 ) THEN
@@ -523,13 +523,13 @@ CONTAINS
DO iz = 1,local_nz
DO iky = 1,local_nky
DO ikx = 1,local_nkx
- Tperp_ = 0._dp
+ Tperp_ = 0._xp
DO ij = 1, local_nj
- j_dp = jarray(ij+ngj/2)
+ j_xp = jarray(ij+ngj/2)
Tperp_ = Tperp_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven)*&
- ((2_dp*j_dp+1)*moments(ia,ip0,ij +ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
- -j_dp*moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
- -(j_dp+1)*moments(ia,ip0,ij+1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
+ ((2_xp*j_xp+1)*moments(ia,ip0,ij +ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
+ -j_xp*moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
+ -(j_xp+1)*moments(ia,ip0,ij+1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
ENDDO
Tper(ia,iky,ikx,iz) = Tperp_
ENDDO
@@ -543,8 +543,8 @@ CONTAINS
SUBROUTINE compute_Tpar
USE time_integration, ONLY : updatetlevel
IMPLICIT NONE
- REAL(dp) :: j_dp
- COMPLEX(dp) :: Tpar_
+ REAL(xp) :: j_xp
+ COMPLEX(xp) :: Tpar_
INTEGER :: ia, iz, iky, ikx, ij
DO ia=1,local_na
IF ( CONTAINSp0 .AND. CONTAINSp0 ) THEN
@@ -552,9 +552,9 @@ CONTAINS
DO iz = 1,local_nz
DO iky = 1,local_nky
DO ikx = 1,local_nkx
- Tpar_ = 0._dp
+ Tpar_ = 0._xp
DO ij = 1, local_nj
- j_dp = REAL(ij-1,dp)
+ j_xp = REAL(ij-1,xp)
Tpar_ = Tpar_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven)*&
(SQRT2 * moments(ia,ip2,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel) &
+ moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
@@ -576,7 +576,7 @@ CONTAINS
CALL compute_Tpar
CALL compute_Tperp
! Temperature
- temp = (Tpar - 2._dp * Tper)/3._dp - dens
+ temp = (Tpar - 2._xp * Tper)/3._xp - dens
END SUBROUTINE compute_fluid_moments
END MODULE processing
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index 4a574d67..142231aa 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -19,15 +19,15 @@ CONTAINS
USE parallel, ONLY: comm0
IMPLICIT NONE
CHARACTER(LEN=50) :: dset_name
- REAL(dp):: time_cp
+ REAL(xp):: time_cp
INTEGER :: cstep_cp, jobnum_cp
INTEGER :: n_
INTEGER :: deltap_cp
INTEGER :: pmax_cp, jmax_cp, n0, Nkx_cp, Nky_cp, Nz_cp, Na_cp, Np_cp, Nj_cp
INTEGER :: ia,ip,ij,iky,ikx,iz, iacp,ipcp,ijcp,iycp,ixcp,izcp, ierr
INTEGER :: ipi,iji,izi
- REAL(dp):: timer_tot_1,timer_tot_2
- COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: moments_cp
+ REAL(xp):: timer_tot_1,timer_tot_2
+ COMPLEX(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: moments_cp
CALL cpu_time(timer_tot_1)
! Checkpoint filename
WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',job2load,'.h5'
@@ -72,7 +72,7 @@ CONTAINS
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments", n_
CALL getarr(fidrst, dset_name, moments_cp(:,:,:,:,:,:))
- moments = 0._dp;
+ moments = 0._xp;
z: DO iz = 1,local_nz
izcp = iz + local_nz_offset
izi = iz + ngz/2
diff --git a/src/solve_EM_fields.F90 b/src/solve_EM_fields.F90
index 73464ae2..01fc26ff 100644
--- a/src/solve_EM_fields.F90
+++ b/src/solve_EM_fields.F90
@@ -7,7 +7,7 @@ SUBROUTINE solve_EM_fields
CALL poisson
- IF(beta .GT. 0._dp) &
+ IF(beta .GT. 0._xp) &
CALL ampere
CONTAINS
@@ -29,8 +29,8 @@ CONTAINS
IMPLICIT NONE
INTEGER :: in, ia, ikx, iky, iz, izi, ini
- COMPLEX(dp) :: fsa_phi, intf_, rhtot ! current flux averaged phi
- COMPLEX(dp), DIMENSION(local_nz) :: rho, integrant ! charge density q_a n_a and aux var
+ COMPLEX(xp) :: fsa_phi, intf_, rhtot ! current flux averaged phi
+ COMPLEX(xp), DIMENSION(local_nz) :: rho, integrant ! charge density q_a n_a and aux var
rhtot = 0
!! Poisson can be solved only for process containng p=0
IF ( SOLVE_POISSON ) THEN
@@ -39,7 +39,7 @@ CONTAINS
!!!!!!!!!!!!!!! Compute particle charge density q_a n_a for each evolved species
DO iz = 1,local_nz
izi = iz+ngz/2
- rho(iz) = 0._dp
+ rho(iz) = 0._xp
DO in = 1,total_nj
ini = in+ngj/2
DO ia = 1,local_na
@@ -55,8 +55,8 @@ CONTAINS
! [qi^2(1-sum_j K_j^2)/tau_i] <phi>_psi = <q_i n_i >_psi
! inv_pol_ion^-1 fsa_phi = simpson(Jacobian rho_i ) * iInt_Jacobian
IF (ADIAB_E) THEN
- fsa_phi = 0._dp
- IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (y-average)
+ fsa_phi = 0._xp
+ IF(kyarray(iky).EQ.0._xp) THEN ! take ky=0 mode (y-average)
! Prepare integrant for z-average
DO iz = 1,local_nz
izi = iz+ngz/2
@@ -78,7 +78,7 @@ CONTAINS
ENDDO y
ENDDO x
! Cancel origin singularity
- IF (contains_kx0 .AND. contains_ky0) phi(iky0,ikx0,:) = 0._dp
+ IF (contains_kx0 .AND. contains_ky0) phi(iky0,ikx0,:) = 0._xp
ENDIF
! Transfer phi to all the others process along p
CALL manual_3D_bcast(phi,local_nky,local_nkx,local_nz+ngz)
@@ -105,7 +105,7 @@ CONTAINS
use species, ONLY: sqrt_tau_o_sigma, q
use model, ONLY: beta
IMPLICIT NONE
- COMPLEX(dp) :: j_a ! current density
+ COMPLEX(xp) :: j_a ! current density
INTEGER :: in, ia, ikx, iky, iz, ini, izi
!! Ampere can be solved only with beta > 0 and for process containng p=1 moments
IF ( SOLVE_AMPERE ) THEN
@@ -114,7 +114,7 @@ CONTAINS
x:DO ikx = 1,local_nkx
y:DO iky = 1,local_nky
!!!!!!!!!!!!!!! compute current density contribution "j_a = q_a u_a" for each species
- j_a = 0._dp
+ j_a = 0._xp
n:DO in=1,total_nj
ini = in+ngj/2
a:DO ia = 1,local_na
@@ -129,7 +129,7 @@ CONTAINS
ENDDO z
ENDIF
! Cancel origin singularity
- IF (contains_kx0 .AND. contains_ky0) psi(iky0,ikx0,:) = 0._dp
+ IF (contains_kx0 .AND. contains_ky0) psi(iky0,ikx0,:) = 0._xp
! Transfer phi to all the others process along p
CALL manual_3D_bcast(psi,local_nky,local_nkx,local_nz+ngz)
END SUBROUTINE ampere
diff --git a/src/species_mod.F90 b/src/species_mod.F90
index 6b40dab8..d00db354 100644
--- a/src/species_mod.F90
+++ b/src/species_mod.F90
@@ -5,31 +5,31 @@ MODULE species
PRIVATE
!! Input parameters
CHARACTER(len=32) :: name_ ! name of the species
- REAL(dp) :: tau_ ! Temperature
- REAL(dp) :: sigma_ ! sqrt mass ratio
- REAL(dp) :: q_ ! Charge
- REAL(dp) :: k_N_ ! density drive (L_ref/L_Ni)
- REAL(dp) :: k_T_ ! temperature drive (L_ref/L_Ti)
+ REAL(xp) :: tau_ ! Temperature
+ REAL(xp) :: sigma_ ! sqrt mass ratio
+ REAL(xp) :: q_ ! Charge
+ REAL(xp) :: k_N_ ! density drive (L_ref/L_Ni)
+ REAL(xp) :: k_T_ ! temperature drive (L_ref/L_Ti)
!! Arrays to store all species features
CHARACTER(len=32),&
ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: name ! name of the species
- REAL(dp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: tau ! Temperature
- REAL(dp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: sigma ! sqrt mass ratio
- REAL(dp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: q ! Charge
- REAL(dp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: k_N ! density drive (L_ref/L_Ni)
- REAL(dp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: k_T ! temperature drive (L_ref/L_Ti)
- REAL(dp), ALLOCATABLE, DIMENSION(:,:),PUBLIC, PROTECTED :: nu_ab ! Collision frequency tensor
+ REAL(xp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: tau ! Temperature
+ REAL(xp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: sigma ! sqrt mass ratio
+ REAL(xp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: q ! Charge
+ REAL(xp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: k_N ! density drive (L_ref/L_Ni)
+ REAL(xp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: k_T ! temperature drive (L_ref/L_Ti)
+ REAL(xp), ALLOCATABLE, DIMENSION(:,:),PUBLIC, PROTECTED :: nu_ab ! Collision frequency tensor
!! Auxiliary variables to store precomputation
- REAL(dp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: tau_q ! factor of the magnetic moment coupling
- REAL(dp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: q_tau !
- REAL(dp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: sqrtTau_q ! factor of parallel moment term
- REAL(dp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: q_sigma_sqrtTau ! factor of parallel phi term
- REAL(dp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: sigma2_tau_o2 ! factor of the Kernel argument
- REAL(dp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: sqrt_sigma2_tau_o2 ! to avoid multiple SQRT eval
- REAL(dp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: q2_tau ! factor of the gammaD sum
- REAL(dp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: q_o_sqrt_tau_sigma ! For psi field terms
- REAL(dp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: sqrt_tau_o_sigma ! For Ampere eq
- REAL(dp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: dpdx ! radial pressure gradient
+ REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: tau_q ! factor of the magnetic moment coupling
+ REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: q_tau !
+ REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: sqrtTau_q ! factor of parallel moment term
+ REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: q_sigma_sqrtTau ! factor of parallel phi term
+ REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: sigma2_tau_o2 ! factor of the Kernel argument
+ REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: sqrt_sigma2_tau_o2 ! to avoid multiple SQRT eval
+ REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: q2_tau ! factor of the gammaD sum
+ REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: q_o_sqrt_tau_sigma ! For psi field terms
+ REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: sqrt_tau_o_sigma ! For Ampere eq
+ REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: xpdx ! radial pressure gradient
!! Accessible routines
PUBLIC :: species_readinputs, species_outputinputs
CONTAINS
@@ -50,11 +50,11 @@ CONTAINS
DO ia = 1,Na
! default parameters
name_ = 'ions'
- tau_ = 1._dp
- sigma_ = 1._dp
- q_ = 1._dp
- k_N_ = 2.22_dp
- k_T_ = 6.96_dp
+ tau_ = 1._xp
+ sigma_ = 1._xp
+ q_ = 1._xp
+ k_N_ = 2.22_xp
+ k_T_ = 6.96_xp
! read input
READ(lu_in,species)
! place values found in the arrays
@@ -69,12 +69,12 @@ CONTAINS
q_tau(ia) = q_/tau_
sqrtTau_q(ia) = sqrt(tau_)/q_
q_sigma_sqrtTau(ia) = q_/sigma_/SQRT(tau_)
- sigma2_tau_o2(ia) = sigma_**2 * tau_/2._dp
- sqrt_sigma2_tau_o2(ia) = SQRT(sigma_**2 * tau_/2._dp)
+ sigma2_tau_o2(ia) = sigma_**2 * tau_/2._xp
+ sqrt_sigma2_tau_o2(ia) = SQRT(sigma_**2 * tau_/2._xp)
q2_tau(ia) = (q_**2)/tau_
q_o_sqrt_tau_sigma(ia) = q_/SQRT(tau_)/sigma_
sqrt_tau_o_sigma(ia) = SQRT(tau_)/sigma_
- dpdx(ia) = 0._dp !not implemented yet
+ xpdx(ia) = 0._xp !not implemented yet
! We remove the adiabatic electron flag if electrons are included
SELECT CASE (name_)
CASE ('electrons','e','electron')
@@ -91,18 +91,18 @@ CONTAINS
! nu_ii = 4 sqrt(pi)/3 T_i^(-3/2) m_i^(-1/2) q^4 n_i0 ln(Lambda)
SELECT CASE (name(ia))
CASE ('electrons','e','electron') ! e-e and e-i collision
- nu_ab(ia,ib) = nu/sigma(ia) * (tau(ia))**(3._dp/2._dp) ! (where already multiplied by 0.532)
+ nu_ab(ia,ib) = nu/sigma(ia) * (tau(ia))**(3._xp/2._xp) ! (where already multiplied by 0.532)
CASE ('ions','ion','i') ! i-e and i-i collision
nu_ab(ia,ib) = nu
CASE DEFAULT
ERROR STOP "!! No collision model for these species interactions"
END SELECT
! I think we can just write
- ! nu_ab(ia,ib) = nu/sigma(ia) * (tau(ia))**(3._dp/2._dp)
+ ! nu_ab(ia,ib) = nu/sigma(ia) * (tau(ia))**(3._xp/2._xp)
ENDDO
ENDDO
ENDIF
- ! nu_e = nu/sigma_e * (tau_e)**(3._dp/2._dp) ! electron-ion collision frequency (where already multiplied by 0.532)
+ ! nu_e = nu/sigma_e * (tau_e)**(3._xp/2._xp) ! electron-ion collision frequency (where already multiplied by 0.532)
! nu_i = nu ! ion-ion collision frequ.
! nu_ee = nu_e ! e-e coll. frequ.
! nu_ie = nu_i ! i-e coll. frequ.
@@ -150,7 +150,7 @@ CONTAINS
ALLOCATE( q2_tau(1:Na))
ALLOCATE(q_o_sqrt_tau_sigma(1:Na))
ALLOCATE( sqrt_tau_o_sigma(1:Na))
- ALLOCATE( dpdx(1:Na))
+ ALLOCATE( xpdx(1:Na))
END SUBROUTINE species_allocate
END MODULE species
diff --git a/src/stepon.F90 b/src/stepon.F90
index e145ecd6..f388a649 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -7,7 +7,7 @@ SUBROUTINE stepon
USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
use mpi, ONLY: MPI_COMM_WORLD
USE time_integration, ONLY: ntimelevel
- USE prec_const, ONLY: dp
+ USE prec_const, ONLY: xp
IMPLICIT NONE
INTEGER :: num_step, ierr
@@ -158,7 +158,7 @@ CONTAINS
END DO
! must be real at origin and top right
moments(:,:,:, iky0,ikx0,:,updatetlevel) = &
- REAL(moments(:,:,:, iky0,ikx0,:,updatetlevel),dp)
+ REAL(moments(:,:,:, iky0,ikx0,:,updatetlevel),xp)
! ENDDO a
! ENDDO p
! ENDDO j
@@ -168,13 +168,13 @@ CONTAINS
phi(iky0,ikx,:) = phi(iky0,total_nkx+2-ikx,:)
END DO
! must be real at origin
- phi(iky0,ikx0,:) = REAL(phi(iky0,ikx0,:),dp)
+ phi(iky0,ikx0,:) = REAL(phi(iky0,ikx0,:),xp)
! Psi
DO ikx=2,total_nkx/2 !symmetry at ky = 0
psi(iky0,ikx,:) = psi(iky0,total_nkx+2-ikx,:)
END DO
! must be real at origin
- psi(iky0,ikx0,:) = REAL(psi(iky0,ikx0,:),dp)
+ psi(iky0,ikx0,:) = REAL(psi(iky0,ikx0,:),xp)
ENDIF
END SUBROUTINE enforce_symmetry
diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90
index 7f937147..e1a2aa51 100644
--- a/src/time_integration_mod.F90
+++ b/src/time_integration_mod.F90
@@ -7,9 +7,9 @@ MODULE time_integration
INTEGER, PUBLIC, PROTECTED :: ntimelevel=4 ! Total number of time levels required by the numerical scheme
INTEGER, PUBLIC, PROTECTED :: updatetlevel ! Current time level to be updated
- real(dp),PUBLIC,PROTECTED,DIMENSION(:,:),ALLOCATABLE :: A_E,A_I
- real(dp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: b_E,b_Es,b_I
- real(dp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: c_E,c_I !Coeff for Expl/Implic time integration in case of time dependent RHS (i.e. dy/dt = f(y,t)) see Baptiste Frei CSE Rapport 06/17
+ real(xp),PUBLIC,PROTECTED,DIMENSION(:,:),ALLOCATABLE :: A_E,A_I
+ real(xp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: b_E,b_Es,b_I
+ real(xp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: c_E,c_I !Coeff for Expl/Implic time integration in case of time dependent RHS (i.e. dy/dt = f(y,t)) see Baptiste Frei CSE Rapport 06/17
character(len=10),PUBLIC,PROTECTED :: numerical_scheme='RK4'
@@ -88,11 +88,11 @@ CONTAINS
CALL allocate_array(b_E,1,nbstep)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 2
- c_E(1) = 0.0_dp
- c_E(2) = 1.0_dp
- b_E(1) = 1._dp/2._dp
- b_E(2) = 1._dp/2._dp
- A_E(2,1) = 1._dp
+ c_E(1) = 0.0_xp
+ c_E(2) = 1.0_xp
+ b_E(1) = 1._xp/2._xp
+ b_E(2) = 1._xp/2._xp
+ A_E(2,1) = 1._xp
END SUBROUTINE RK2
SUBROUTINE SSPx_RK2
@@ -103,21 +103,21 @@ CONTAINS
USE prec_const
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep = 2
- REAL(dp) :: alpha, beta
- alpha = 1._dp/SQRT(2._dp)
- beta = SQRT(2._dp) - 1._dp
+ REAL(xp) :: alpha, beta
+ alpha = 1._xp/SQRT(2._xp)
+ beta = SQRT(2._xp) - 1._xp
CALL allocate_array(c_E,1,nbstep)
CALL allocate_array(b_E,1,nbstep)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 2
- c_E(1) = 0.0_dp
- c_E(2) = 1.0_dp/2.0_dp
- b_E(1) = alpha*beta/2._dp
- b_E(2) = alpha/2._dp
+ c_E(1) = 0.0_xp
+ c_E(2) = 1.0_xp/2.0_xp
+ b_E(1) = alpha*beta/2._xp
+ b_E(2) = alpha/2._xp
A_E(2,1) = alpha
- ! b_E(1) = 1._dp
- ! b_E(2) = 1._dp/SQRT(2._dp)
- ! A_E(2,1) = 1._dp/SQRT(2._dp)
+ ! b_E(1) = 1._xp
+ ! b_E(2) = 1._xp/SQRT(2._xp)
+ ! A_E(2,1) = 1._xp/SQRT(2._xp)
END SUBROUTINE SSPx_RK2
!!! third order time schemes
@@ -131,15 +131,15 @@ CONTAINS
CALL allocate_array(b_E,1,nbstep)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 3
- c_E(1) = 0.0_dp
- c_E(2) = 1.0_dp/2.0_dp
- c_E(3) = 1.0_dp
- b_E(1) = 1._dp/6._dp
- b_E(2) = 2._dp/3._dp
- b_E(3) = 1._dp/6._dp
- A_E(2,1) = 1.0_dp/2.0_dp
- A_E(3,1) = -1._dp
- A_E(3,2) = 2._dp
+ c_E(1) = 0.0_xp
+ c_E(2) = 1.0_xp/2.0_xp
+ c_E(3) = 1.0_xp
+ b_E(1) = 1._xp/6._xp
+ b_E(2) = 2._xp/3._xp
+ b_E(3) = 1._xp/6._xp
+ A_E(2,1) = 1.0_xp/2.0_xp
+ A_E(3,1) = -1._xp
+ A_E(3,2) = 2._xp
END SUBROUTINE RK3
SUBROUTINE SSPx_RK3
@@ -150,24 +150,24 @@ CONTAINS
USE prec_const
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep = 3
- REAL(dp) :: a1, a2, a3, w1, w2, w3
- a1 = (1._dp/6._dp)**(1._dp/3._dp)! (1/6)^(1/3)
- ! a1 = 0.5503212081491044571635029569733887910843_dp ! (1/6)^(1/3)
+ REAL(xp) :: a1, a2, a3, w1, w2, w3
+ a1 = (1._xp/6._xp)**(1._xp/3._xp)! (1/6)^(1/3)
+ ! a1 = 0.5503212081491044571635029569733887910843_xp ! (1/6)^(1/3)
a2 = a1
a3 = a1
- w1 = 0.5_dp*(-1._dp + SQRT( 9._dp - 2._dp * 6._dp**(2._dp/3._dp))) ! (-1 + sqrt(9-2*6^(2/3)))/2
- ! w1 = 0.2739744023885328783052273138309828937054_dp ! (sqrt(9-2*6^(2/3))-1)/2
- w2 = 0.5_dp*(-1._dp + 6._dp**(2._dp/3._dp) - SQRT(9._dp - 2._dp * 6._dp**(2._dp/3._dp))) ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2
- ! w2 = 0.3769892220587804931852815570891834795475_dp ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2
- w3 = 1._dp/a1 - w2 * (1._dp + w1)
- ! w3 = 1.3368459739528868457369981115334667265415_dp
+ w1 = 0.5_xp*(-1._xp + SQRT( 9._xp - 2._xp * 6._xp**(2._xp/3._xp))) ! (-1 + sqrt(9-2*6^(2/3)))/2
+ ! w1 = 0.2739744023885328783052273138309828937054_xp ! (sqrt(9-2*6^(2/3))-1)/2
+ w2 = 0.5_xp*(-1._xp + 6._xp**(2._xp/3._xp) - SQRT(9._xp - 2._xp * 6._xp**(2._xp/3._xp))) ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2
+ ! w2 = 0.3769892220587804931852815570891834795475_xp ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2
+ w3 = 1._xp/a1 - w2 * (1._xp + w1)
+ ! w3 = 1.3368459739528868457369981115334667265415_xp
CALL allocate_array(c_E,1,nbstep)
CALL allocate_array(b_E,1,nbstep)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 3
- c_E(1) = 0.0_dp
- c_E(2) = 1.0_dp/2.0_dp
- c_E(3) = 1.0_dp/2.0_dp
+ c_E(1) = 0.0_xp
+ c_E(2) = 1.0_xp/2.0_xp
+ c_E(3) = 1.0_xp/2.0_xp
b_E(1) = a1 * (w1*w2 + w3)
b_E(2) = a2 * w2
b_E(3) = a3
@@ -187,15 +187,15 @@ CONTAINS
CALL allocate_array(b_E,1,nbstep)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 3
- c_E(1) = 0._dp
- c_E(2) = 0.711664700366941_dp
- c_E(3) = 0.994611536833690_dp
- b_E(1) = 0.398930808264688_dp
- b_E(2) = 0.345755244189623_dp
- b_E(3) = 0.255313947545689_dp
- A_E(2,1) = 0.711664700366941_dp
- A_E(3,1) = 0.077338168947683_dp
- A_E(3,2) = 0.917273367886007_dp
+ c_E(1) = 0._xp
+ c_E(2) = 0.711664700366941_xp
+ c_E(3) = 0.994611536833690_xp
+ b_E(1) = 0.398930808264688_xp
+ b_E(2) = 0.345755244189623_xp
+ b_E(3) = 0.255313947545689_xp
+ A_E(2,1) = 0.711664700366941_xp
+ A_E(3,1) = 0.077338168947683_xp
+ A_E(3,2) = 0.917273367886007_xp
END SUBROUTINE IMEX_SSP2
SUBROUTINE ARK2
@@ -209,15 +209,15 @@ CONTAINS
CALL allocate_array(b_E,1,nbstep)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 3
- c_E(1) = 0._dp
- c_E(2) = 2._dp*(1._dp - 1._dp/SQRT2)
- c_E(3) = 1._dp
- b_E(1) = 1._dp/(2._dp*SQRT2)
- b_E(2) = 1._dp/(2._dp*SQRT2)
- b_E(3) = 1._dp - 1._dp/SQRT2
- A_E(2,1) = 2._dp*(1._dp - 1._dp/SQRT2)
- A_E(3,1) = 1._dp - (3._dp + 2._dp*SQRT2)/6._dp
- A_E(3,2) = (3._dp + 2._dp*SQRT2)/6._dp
+ c_E(1) = 0._xp
+ c_E(2) = 2._xp*(1._xp - 1._xp/SQRT2)
+ c_E(3) = 1._xp
+ b_E(1) = 1._xp/(2._xp*SQRT2)
+ b_E(2) = 1._xp/(2._xp*SQRT2)
+ b_E(3) = 1._xp - 1._xp/SQRT2
+ A_E(2,1) = 2._xp*(1._xp - 1._xp/SQRT2)
+ A_E(3,1) = 1._xp - (3._xp + 2._xp*SQRT2)/6._xp
+ A_E(3,2) = (3._xp + 2._xp*SQRT2)/6._xp
END SUBROUTINE ARK2
SUBROUTINE SSP_RK3
@@ -230,15 +230,15 @@ CONTAINS
CALL allocate_array(b_E,1,nbstep)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 3
- c_E(1) = 0.0_dp
- c_E(2) = 1.0_dp
- c_E(3) = 1.0_dp/2.0_dp
- b_E(1) = 1._dp/6._dp
- b_E(2) = 1._dp/6._dp
- b_E(3) = 2._dp/3._dp
- A_E(2,1) = 1._dp
- A_E(3,1) = 1._dp/4._dp
- A_E(3,2) = 1._dp/4._dp
+ c_E(1) = 0.0_xp
+ c_E(2) = 1.0_xp
+ c_E(3) = 1.0_xp/2.0_xp
+ b_E(1) = 1._xp/6._xp
+ b_E(2) = 1._xp/6._xp
+ b_E(3) = 2._xp/3._xp
+ A_E(2,1) = 1._xp
+ A_E(3,1) = 1._xp/4._xp
+ A_E(3,2) = 1._xp/4._xp
END SUBROUTINE SSP_RK3
!!! fourth order time schemes
@@ -252,17 +252,17 @@ CONTAINS
CALL allocate_array(b_E,1,nbstep)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 4
- c_E(1) = 0.0_dp
- c_E(2) = 1.0_dp/2.0_dp
- c_E(3) = 1.0_dp/2.0_dp
- c_E(4) = 1.0_dp
- b_E(1) = 1.0_dp/6.0_dp
- b_E(2) = 1.0_dp/3.0_dp
- b_E(3) = 1.0_dp/3.0_dp
- b_E(4) = 1.0_dp/6.0_dp
- A_E(2,1) = 1.0_dp/2.0_dp
- A_E(3,2) = 1.0_dp/2.0_dp
- A_E(4,3) = 1.0_dp
+ c_E(1) = 0.0_xp
+ c_E(2) = 1.0_xp/2.0_xp
+ c_E(3) = 1.0_xp/2.0_xp
+ c_E(4) = 1.0_xp
+ b_E(1) = 1.0_xp/6.0_xp
+ b_E(2) = 1.0_xp/3.0_xp
+ b_E(3) = 1.0_xp/3.0_xp
+ b_E(4) = 1.0_xp/6.0_xp
+ A_E(2,1) = 1.0_xp/2.0_xp
+ A_E(3,2) = 1.0_xp/2.0_xp
+ A_E(4,3) = 1.0_xp
END SUBROUTINE RK4
!!! fifth order time schemes
@@ -277,40 +277,40 @@ CONTAINS
CALL allocate_array(b_E,1,nbstep)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 7
- c_E(1) = 0._dp
- c_E(2) = 1.0_dp/5.0_dp
- c_E(3) = 3.0_dp /10.0_dp
- c_E(4) = 4.0_dp/5.0_dp
- c_E(5) = 8.0_dp/9.0_dp
- c_E(6) = 1.0_dp
- c_E(7) = 1.0_dp
- A_E(2,1) = 1.0_dp/5.0_dp
- A_E(3,1) = 3.0_dp/40.0_dp
- A_E(3,2) = 9.0_dp/40.0_dp
- A_E(4,1) = 44.0_dp/45.0_dp
- A_E(4,2) = -56.0_dp/15.0_dp
- A_E(4,3) = 32.0_dp/9.0_dp
- A_E(5,1 ) = 19372.0_dp/6561.0_dp
- A_E(5,2) = -25360.0_dp/2187.0_dp
- A_E(5,3) = 64448.0_dp/6561.0_dp
- A_E(5,4) = -212.0_dp/729.0_dp
- A_E(6,1) = 9017.0_dp/3168.0_dp
- A_E(6,2)= -355.0_dp/33.0_dp
- A_E(6,3) = 46732.0_dp/5247.0_dp
- A_E(6,4) = 49.0_dp/176.0_dp
- A_E(6,5) = -5103.0_dp/18656.0_dp
- A_E(7,1) = 35.0_dp/384.0_dp
- A_E(7,3) = 500.0_dp/1113.0_dp
- A_E(7,4) = 125.0_dp/192.0_dp
- A_E(7,5) = -2187.0_dp/6784.0_dp
- A_E(7,6) = 11.0_dp/84.0_dp
- b_E(1) = 35.0_dp/384.0_dp
- b_E(2) = 0._dp
- b_E(3) = 500.0_dp/1113.0_dp
- b_E(4) = 125.0_dp/192.0_dp
- b_E(5) = -2187.0_dp/6784.0_dp
- b_E(6) = 11.0_dp/84.0_dp
- b_E(7) = 0._dp
+ c_E(1) = 0._xp
+ c_E(2) = 1.0_xp/5.0_xp
+ c_E(3) = 3.0_xp /10.0_xp
+ c_E(4) = 4.0_xp/5.0_xp
+ c_E(5) = 8.0_xp/9.0_xp
+ c_E(6) = 1.0_xp
+ c_E(7) = 1.0_xp
+ A_E(2,1) = 1.0_xp/5.0_xp
+ A_E(3,1) = 3.0_xp/40.0_xp
+ A_E(3,2) = 9.0_xp/40.0_xp
+ A_E(4,1) = 44.0_xp/45.0_xp
+ A_E(4,2) = -56.0_xp/15.0_xp
+ A_E(4,3) = 32.0_xp/9.0_xp
+ A_E(5,1 ) = 19372.0_xp/6561.0_xp
+ A_E(5,2) = -25360.0_xp/2187.0_xp
+ A_E(5,3) = 64448.0_xp/6561.0_xp
+ A_E(5,4) = -212.0_xp/729.0_xp
+ A_E(6,1) = 9017.0_xp/3168.0_xp
+ A_E(6,2)= -355.0_xp/33.0_xp
+ A_E(6,3) = 46732.0_xp/5247.0_xp
+ A_E(6,4) = 49.0_xp/176.0_xp
+ A_E(6,5) = -5103.0_xp/18656.0_xp
+ A_E(7,1) = 35.0_xp/384.0_xp
+ A_E(7,3) = 500.0_xp/1113.0_xp
+ A_E(7,4) = 125.0_xp/192.0_xp
+ A_E(7,5) = -2187.0_xp/6784.0_xp
+ A_E(7,6) = 11.0_xp/84.0_xp
+ b_E(1) = 35.0_xp/384.0_xp
+ b_E(2) = 0._xp
+ b_E(3) = 500.0_xp/1113.0_xp
+ b_E(4) = 125.0_xp/192.0_xp
+ b_E(5) = -2187.0_xp/6784.0_xp
+ b_E(6) = 11.0_xp/84.0_xp
+ b_E(7) = 0._xp
END SUBROUTINE DOPRI5
END MODULE time_integration
diff --git a/src/utility_mod.F90 b/src/utility_mod.F90
index 601b6b7f..b8544848 100644
--- a/src/utility_mod.F90
+++ b/src/utility_mod.F90
@@ -6,7 +6,7 @@ CONTAINS
FUNCTION is_nan(x,str) RESULT(isn)
USE basic, ONLY: cstep
USE time_integration, ONLY: updatetlevel
- USE prec_const, ONLY: dp, stdout
+ USE prec_const, ONLY: xp, stdout
IMPLICIT NONE
real, INTENT(IN) :: x
@@ -25,7 +25,7 @@ CONTAINS
FUNCTION is_inf(x,str) RESULT(isi)
- USE prec_const, ONLY: dp, stdout
+ USE prec_const, ONLY: xp, stdout
IMPLICIT NONE
real, INTENT(IN) :: x
@@ -44,17 +44,17 @@ CONTAINS
END FUNCTION is_inf
! FUNCTION checkfield(n1,n2,n3,field,str) RESULT(mlend)
- ! use prec_const, ONLY: dp
+ ! use prec_const, ONLY: xp
! IMPLICIT NONE
! !! BUG found (or feature?)
! ! if one put the commented first line (commented) instead of the second one,
! ! no error will be risen by the compiler even if the rank of the array is not matching (should be 3D!)
- ! ! COMPLEX(dp), DIMENSION(ikys:ikye,ikxs:ikxe), INTENT(IN) :: field
+ ! ! COMPLEX(xp), DIMENSION(ikys:ikye,ikxs:ikxe), INTENT(IN) :: field
! INTEGER, INTENT(in) :: n1,n2,n3
- ! COMPLEX(dp), DIMENSION(n1,n2,n3), INTENT(IN) :: field
+ ! COMPLEX(xp), DIMENSION(n1,n2,n3), INTENT(IN) :: field
! CHARACTER(LEN=*), INTENT(IN) :: str
! LOGICAL :: mlend
- ! COMPLEX(dp) :: sumfield
+ ! COMPLEX(xp) :: sumfield
! sumfield=SUM(field)
@@ -63,9 +63,9 @@ CONTAINS
! END FUNCTION checkfield
! FUNCTION checkelem(elem,str) RESULT(mlend)
- ! use prec_const, ONLY: dp
+ ! use prec_const, ONLY: xp
! IMPLICIT NONE
- ! COMPLEX(dp), INTENT(IN) :: elem
+ ! COMPLEX(xp), INTENT(IN) :: elem
! CHARACTER(LEN=*), INTENT(IN) :: str
! LOGICAL :: mlend
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index fdd59170..34d685b3 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -5,8 +5,13 @@ PARTITION = '/misc/gyacomo23_outputs/';
%% CBC
% resdir = 'paper_2_GYAC23/CBC/7x4x192x96x32_nu_0.05_muxy_0.5_muz_0.2';
-resdir = 'paper_2_GYAC23/CBC/7x4x192x96x32_nu_0.05_muxy_1.0_muz_1.0;'
+% resdir = 'paper_2_GYAC23/CBC/7x4x192x96x32_nu_0.05_muxy_1.0_muz_1.0';
+% resdir = 'paper_2_GYAC23/CBC/7x4x192x96x32_nu_0.05_muxy_1.0_muz_2.0';
+% resdir = 'paper_2_GYAC23/CBC/Full_NL_7x4x192x96x32_nu_0.05_muxy_1.0_muz_2.0';
+%% tests
+resdir = 'paper_2_GYAC23/precision_study/5x3x128x64x24';
+% resdir = 'paper_2_GYAC23/precision_study/5x3x128x64x24_Lx_180';
%%
J0 = 00; J1 = 10;
@@ -28,6 +33,7 @@ options.INTERP = 0;
options.RESOLUTION = 256;
fig = plot_radial_transport_and_spacetime(data,options);
+if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
options.INTERP = 0;
@@ -40,7 +46,7 @@ options.NAME = '\phi';
% options.NAME = 'Q_x';
% options.NAME = 'n_i';
% options.NAME = 'n_i-n_e';
-options.PLAN = 'xz';
+options.PLAN = 'xy';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
@@ -51,3 +57,4 @@ data.EPS = 0.1;
data.a = data.EPS * 2000;
options.RESOLUTION = 256;
create_film(data,options,'.gif')
+end
--
GitLab