From 2532ff3a96ce0b7a5891d05224d64fbcc41c04d2 Mon Sep 17 00:00:00 2001 From: Andreas Zach Date: Sat, 10 Jan 2026 00:42:22 +0100 Subject: [PATCH 1/3] Make all reals real(dp) where dp => real64 --- src/attenuation_factor.f90 | 20 ++-- src/collis_nbi.f90 | 96 ++++++++-------- src/datatypes.f90 | 113 +++++++++---------- src/diag/diag_atten_map.f90 | 38 +++---- src/diag/diag_bounce_debug.f90 | 76 ++++++------- src/diag/diag_bounce_nonlin.f90 | 36 +++--- src/diag/diag_contrib_map.f90 | 80 +++++++------- src/do_magfie_standalone.f90 | 188 ++++++++++++++++---------------- src/driftorbit.f90 | 21 ++-- src/freq.f90 | 145 ++++++++++++------------ src/magfie.f90 | 29 ++--- src/neort.f90 | 108 +++++++++--------- src/nonlin.f90 | 25 +++-- src/orbit.f90 | 182 ++++++++++++++++--------------- src/profiles.f90 | 38 +++---- src/resonance.f90 | 51 ++++----- src/transport.f90 | 85 ++++++++------- src/util.f90 | 42 +++---- test/test_bounce.f90 | 45 ++++---- test/test_frequencies.f90 | 19 ++-- test/test_misc.f90 | 67 ++++++------ test/test_neort_lib.f90 | 18 +-- test/test_omega_prime.f90 | 131 +++++++++++----------- test/test_parallel.f90 | 8 +- test/test_reslines.f90 | 21 ++-- test/test_torque.f90 | 55 +++++----- 26 files changed, 883 insertions(+), 854 deletions(-) diff --git a/src/attenuation_factor.f90 b/src/attenuation_factor.f90 index d02c7c91..1e258805 100644 --- a/src/attenuation_factor.f90 +++ b/src/attenuation_factor.f90 @@ -1,10 +1,12 @@ module thetadata_mod + use iso_fortran_env, only: dp => real64 + implicit none logical :: attenuation_data_initialized = .false. integer :: npoiarg - double precision :: argmin, argstepinv - double precision, dimension(:), allocatable :: argvals, thetavals + real(dp) :: argmin, argstepinv + real(dp), dimension(:), allocatable :: argvals, thetavals ! All variables are shared (read-only after init_attenuation_data) @@ -15,7 +17,7 @@ subroutine init_attenuation_data() ! MUST be called BEFORE any parallel region integer, parameter :: iunitarg = 741 integer :: i - double precision :: arglog, theta_val + real(dp) :: arglog, theta_val if (attenuation_data_initialized) return @@ -40,7 +42,7 @@ subroutine init_attenuation_data() close (iunitarg) argmin = argvals(1) - argstepinv = 1.d0/(argvals(2) - argvals(1)) + argstepinv = 1.0_dp/(argvals(2) - argvals(1)) attenuation_data_initialized = .true. end subroutine init_attenuation_data @@ -48,17 +50,17 @@ end subroutine init_attenuation_data end module thetadata_mod subroutine attenuation_factor(D, Theta) - + use iso_fortran_env, only: dp => real64 use thetadata_mod, only: npoiarg, argmin, argstepinv, argvals, thetavals use polylag_3, only: mp, indef, plag1d implicit none - double precision, intent(in) :: D - double precision, intent(out) :: Theta - double precision :: arglog, fun, der + real(dp), intent(in) :: D + real(dp), intent(out) :: Theta + real(dp) :: arglog, fun, der integer, dimension(mp) :: indu - double precision, dimension(mp) :: xp, fp + real(dp), dimension(mp) :: xp, fp arglog = log(D) diff --git a/src/collis_nbi.f90 b/src/collis_nbi.f90 index 4524fee6..be22e44d 100644 --- a/src/collis_nbi.f90 +++ b/src/collis_nbi.f90 @@ -1,8 +1,12 @@ module collis_alp + use iso_fortran_env, only: dp => real64 + + implicit none + integer, parameter :: nsorts = 3, ns = 10000 logical :: swcoll = .false. - double precision, dimension(nsorts) :: efcolf, velrat, enrat + real(dp), dimension(nsorts) :: efcolf, velrat, enrat !$omp threadprivate (swcoll, efcolf, velrat, enrat) @@ -29,24 +33,24 @@ subroutine coleff(p, dpp, dhh, fpeff) implicit none integer i - double precision, intent(in) :: p - double precision, intent(out) :: dpp, dhh, fpeff - double precision :: plim, xbeta, dp, dh, dpd + real(dp), intent(in) :: p + real(dp), intent(out) :: dpp, dhh, fpeff + real(dp) :: plim, xbeta, d_p, dh, dpd - plim = max(p, 1.d-8) + plim = max(p, 1.0e-8_dp) - dpp = 0.0d0 - dhh = 0.0d0 - fpeff = 0.0d0 + dpp = 0.0_dp + dhh = 0.0_dp + fpeff = 0.0_dp do i = 1, nsorts xbeta = p*velrat(i) - call onseff(xbeta, dp, dh, dpd) + call onseff(xbeta, d_p, dh, dpd) - dpp = dpp + dp*efcolf(i) + dpp = dpp + d_p*efcolf(i) dhh = dhh + dh*efcolf(i) - fpeff = fpeff + (dpd/plim - 2.0d0*dp*p*enrat(i))*efcolf(i) + fpeff = fpeff + (dpd/plim - 2.0_dp*d_p*p*enrat(i))*efcolf(i) end do dhh = dhh/plim**2 @@ -54,37 +58,37 @@ subroutine coleff(p, dpp, dhh, fpeff) return end subroutine coleff - subroutine onseff(v, dp, dh, dpd) - ! dp - dimensionless dpp + subroutine onseff(v, d_p, dh, dpd) + ! d_p - dimensionless dpp ! dh - dhh*p^2 (p - dmls) - ! dpd - (1/p)(d/dp)p^2*dp (p - dmls) + ! dpd - (1/p)(d/dp)p^2*d_p (p - dmls) implicit none ! square root of pi - double precision, parameter :: sqp = 1.7724538d0 + real(dp), parameter :: sqp = 1.7724538_dp ! cons=4./(3.*sqrt(pi)) - double precision, parameter :: cons = .75225278d0 - double precision, intent(in) :: v - double precision, intent(out) :: dp, dh, dpd - double precision :: v2, v3, ex, er + real(dp), parameter :: cons = 0.75225278_dp + real(dp), intent(in) :: v + real(dp), intent(out) :: d_p, dh, dpd + real(dp) :: v2, v3, ex, er v2 = v**2 v3 = v2*v - if (v .lt. 0.01d0) then - dp = cons*(1.d0 - 0.6d0*v2) - dh = cons*(1.d0 - 0.2d0*v2) - dpd = 2.d0*cons*(1.d0 - 1.2d0*v2) - elseif (v .gt. 6.d0) then - dp = 1.d0/v3 - dh = (1.d0 - 0.5d0/v2)/v - dpd = -1.d0/v3 + if (v < 0.01_dp) then + d_p = cons*(1.0_dp - 0.6_dp*v2) + dh = cons*(1.0_dp - 0.2_dp*v2) + dpd = 2.0_dp*cons*(1.0_dp - 1.2_dp*v2) + elseif (v > 6.0_dp) then + d_p = 1.0_dp/v3 + dh = (1.0_dp - 0.5_dp/v2)/v + dpd = -1.0_dp/v3 else ex = exp(-v2)/sqp er = erf(v) - dp = er/v3 - 2.d0*ex/v2 - dh = er*(1.d0 - 0.5d0/v2)/v + ex/v2 - dpd = 4.d0*ex - dp + d_p = er/v3 - 2.0_dp*ex/v2 + dh = er*(1.0_dp - 0.5_dp/v2)/v + ex/v2 + dpd = 4.0_dp*ex - d_p end if return @@ -124,36 +128,36 @@ subroutine loacol_nbi(amb, am1, am2, Zb, Z1, Z2, densi1, densi2, tempi1, tempi2, implicit none - double precision :: amb, am1, am2, Zb, Z1, Z2, densi1, densi2, tempi1, tempi2, tempe, ebeam, dense - double precision :: v0, dchichi, slowrate, dchichi_norm, slowrate_norm, vti1, vti2, vte - double precision :: pi, pmass, emass, e, ev, alame, frecol_base, alami1, alami2 + real(dp) :: amb, am1, am2, Zb, Z1, Z2, densi1, densi2, tempi1, tempi2, tempe, ebeam, dense + real(dp) :: v0, dchichi, slowrate, dchichi_norm, slowrate_norm, vti1, vti2, vte + real(dp) :: pi, pmass, emass, e, ev, alame, frecol_base, alami1, alami2 - pi = 3.14159265358979d0 - pmass = 1.6726d-24 - emass = 9.1094d-28 - e = 4.8032d-10 - ev = 1.6022d-12 + pi = 3.14159265358979_dp + pmass = 1.6726e-24_dp + emass = 9.1094e-28_dp + e = 4.8032e-10_dp + ev = 1.6022e-12_dp enrat(1) = ebeam/tempi1 enrat(2) = ebeam/tempi2 enrat(3) = ebeam/tempe - v0 = sqrt(2.d0*ebeam*ev/(amb*pmass)) - vti1 = sqrt(2.d0*tempi1*ev/(pmass*am1)) - vti2 = sqrt(2.d0*tempi2*ev/(pmass*am2)) - vte = sqrt(2.d0*tempe*ev/emass) + v0 = sqrt(2.0_dp*ebeam*ev/(amb*pmass)) + vti1 = sqrt(2.0_dp*tempi1*ev/(pmass*am1)) + vti2 = sqrt(2.0_dp*tempi2*ev/(pmass*am2)) + vte = sqrt(2.0_dp*tempe*ev/emass) velrat(1) = v0/vti1 velrat(2) = v0/vti2 velrat(3) = v0/vte dense = densi1*Z1 + densi2*Z2 - alami1 = 23.d0 - log(max(epsilon(1.d0), & + alami1 = 23.0_dp - log(max(epsilon(1.0_dp), & sqrt(densi1*Z1**2/tempi1)*Zb*Z1*(amb + am1)/(amb*tempi1 + am1*ebeam))) - alami2 = 23.d0 - log(max(epsilon(1.d0), & + alami2 = 23.0_dp - log(max(epsilon(1.0_dp), & sqrt(densi2*Z2**2/tempi2)*Zb*Z2*(amb + am2)/(amb*tempi2 + am2*ebeam))) - alame = 24.d0 - log(sqrt(dense)/tempe) - frecol_base = 2.d0*pi*dense*e**4*Zb**2/((amb*pmass)**2*v0**3) ! usual + alame = 24.0_dp - log(sqrt(dense)/tempe) + frecol_base = 2.0_dp*pi*dense*e**4*Zb**2/((amb*pmass)**2*v0**3) ! usual frecol_base = frecol_base/v0 ! normalized efcolf(1) = frecol_base*Z1**2*alami1*densi1/dense diff --git a/src/datatypes.f90 b/src/datatypes.f90 index 755a99b5..59c0080b 100644 --- a/src/datatypes.f90 +++ b/src/datatypes.f90 @@ -1,49 +1,50 @@ module neort_datatypes + use iso_fortran_env, only: dp => real64 implicit none type :: magfie_params_t - real(8) :: s ! Normalised toroidal flux - real(8) :: R0 ! Major radius [cgs] - real(8) :: a ! Minor radius [cgs] - real(8) :: eps ! Inverse aspect ratio (a/R0) - real(8) :: psi_pr ! Radial derivative of poloidal flux - real(8) :: B0 ! Reference magnetic field magnitude [cgs] - real(8) :: Bthcov ! Covariant poloidal field component - real(8) :: Bphcov ! Covariant toroidal field component - real(8) :: dBthcovds ! Derivative of Bthcov with respect to s - real(8) :: dBphcovds ! Derivative of Bphcov with respect to s - real(8) :: q ! Safety factor - real(8) :: iota ! Rotational transform (1/q) - real(8) :: dVds ! Derivative of flux-surface volume with respect to s - real(8) :: M_t ! Toroidal Mach number - real(8) :: Om_tE ! Toroidal rotation frequency from electric field - real(8) :: Om_tBref ! Reference toroidal rotation frequency - real(8) :: vth ! Thermal velocity [cgs] - real(8) :: T_in_eV ! Temperature [eV] - real(8) :: m0 ! Poloidal mode number of perturbation - real(8) :: n0 ! Toroidal mode number of perturbation - real(8) :: Dp ! Thermodynamic force parameter - real(8) :: Drp_over_Dp ! Ratio of rotation force to thermodynamic force - real(8) :: etatp ! Normalized temperature gradient - real(8) :: etadt ! Normalized density gradient + real(dp) :: s ! Normalised toroidal flux + real(dp) :: R0 ! Major radius [cgs] + real(dp) :: a ! Minor radius [cgs] + real(dp) :: eps ! Inverse aspect ratio (a/R0) + real(dp) :: psi_pr ! Radial derivative of poloidal flux + real(dp) :: B0 ! Reference magnetic field magnitude [cgs] + real(dp) :: Bthcov ! Covariant poloidal field component + real(dp) :: Bphcov ! Covariant toroidal field component + real(dp) :: dBthcovds ! Derivative of Bthcov with respect to s + real(dp) :: dBphcovds ! Derivative of Bphcov with respect to s + real(dp) :: q ! Safety factor + real(dp) :: iota ! Rotational transform (1/q) + real(dp) :: dVds ! Derivative of flux-surface volume with respect to s + real(dp) :: M_t ! Toroidal Mach number + real(dp) :: Om_tE ! Toroidal rotation frequency from electric field + real(dp) :: Om_tBref ! Reference toroidal rotation frequency + real(dp) :: vth ! Thermal velocity [cgs] + real(dp) :: T_in_eV ! Temperature [eV] + real(dp) :: m0 ! Poloidal mode number of perturbation + real(dp) :: n0 ! Toroidal mode number of perturbation + real(dp) :: Dp ! Thermodynamic force parameter + real(dp) :: Drp_over_Dp ! Ratio of rotation force to thermodynamic force + real(dp) :: etatp ! Normalized temperature gradient + real(dp) :: etadt ! Normalized density gradient logical :: pertfile ! Flag: use perturbation from file (not analytic) logical :: nonlin ! Flag: enable nonlinear physics calculations - real(8) :: dpp ! Nonlinear diagnostic parameter - real(8) :: dhh ! Nonlinear diagnostic parameter - real(8) :: fpeff ! Effective passing fraction parameter + real(dp) :: dpp ! Nonlinear diagnostic parameter + real(dp) :: dhh ! Nonlinear diagnostic parameter + real(dp) :: fpeff ! Effective passing fraction parameter end type magfie_params_t type :: magfie_tensors_t - real(8) :: theta ! Boozer poloidal angle [rad] - real(8) :: bmod ! Magnetic field magnitude |B| [cgs] - real(8) :: sqrtg ! Jacobian determinant sqrt(g) - real(8) :: hder(3) ! Contravariant basis vectors - real(8) :: hcovar(3) ! Covariant basis vectors - real(8) :: hctrvr(3) ! Contravariant components of field-line curvature - real(8) :: hcurl(3) ! Components of curl of basis vectors - complex(8) :: bn ! Perturbation amplitude normalised by |B| - complex(8) :: eps_exp ! Reference analytic perturbation epsmn*exp(i*m0*theta) + real(dp) :: theta ! Boozer poloidal angle [rad] + real(dp) :: bmod ! Magnetic field magnitude |B| [cgs] + real(dp) :: sqrtg ! Jacobian determinant sqrt(g) + real(dp) :: hder(3) ! Contravariant basis vectors + real(dp) :: hcovar(3) ! Covariant basis vectors + real(dp) :: hctrvr(3) ! Contravariant components of field-line curvature + real(dp) :: hcurl(3) ! Components of curl of basis vectors + complex(dp) :: bn ! Perturbation amplitude normalised by |B| + complex(dp) :: eps_exp ! Reference analytic perturbation epsmn*exp(i*m0*theta) end type magfie_tensors_t type :: magfie_data_t @@ -54,37 +55,37 @@ module neort_datatypes end type magfie_data_t type :: transport_summary_t - real(8) :: M_t ! Toroidal Mach number - real(8) :: Dco(2) ! Transport coefficients (D11, D12) for co-passing particles - real(8) :: Dctr(2) ! Transport coefficients (D11, D12) for counter-passing particles - real(8) :: Dt(2) ! Transport coefficients (D11, D12) for trapped particles + real(dp) :: M_t ! Toroidal Mach number + real(dp) :: Dco(2) ! Transport coefficients (D11, D12) for co-passing particles + real(dp) :: Dctr(2) ! Transport coefficients (D11, D12) for counter-passing particles + real(dp) :: Dt(2) ! Transport coefficients (D11, D12) for trapped particles ! D11: particle flux coefficient, D12: momentum flux coefficient end type transport_summary_t type :: torque_summary_t logical :: has_torque ! Flag: torque calculation enabled - real(8) :: s ! Normalised toroidal flux - real(8) :: dVds ! Derivative of flux-surface volume with respect to s - real(8) :: M_t ! Toroidal Mach number - real(8) :: Tco ! Torque density from co-passing particles - real(8) :: Tctr ! Torque density from counter-passing particles - real(8) :: Tt ! Torque density from trapped particles + real(dp) :: s ! Normalised toroidal flux + real(dp) :: dVds ! Derivative of flux-surface volume with respect to s + real(dp) :: M_t ! Toroidal Mach number + real(dp) :: Tco ! Torque density from co-passing particles + real(dp) :: Tctr ! Torque density from counter-passing particles + real(dp) :: Tt ! Torque density from trapped particles end type torque_summary_t type :: transport_harmonic_t integer :: mth ! Poloidal mode number associated with resonance - real(8) :: Dresco(2) + real(dp) :: Dresco(2) ! Transport coefficients (D11, D12) for co-passing particles, per harmonic - real(8) :: Dresctr(2) + real(dp) :: Dresctr(2) ! Transport coefficients (D11, D12) for counter-passing particles, per harmonic - real(8) :: Drest(2) ! Transport coefficients (D11, D12) for trapped particles, per harmonic - real(8) :: Tresco ! Torque from co-passing particles, per harmonic - real(8) :: Tresctr ! Torque from counter-passing particles, per harmonic - real(8) :: Trest ! Torque from trapped particles, per harmonic - real(8) :: vminp_over_vth ! Lower bound of passing velocity grid normalised - real(8) :: vmaxp_over_vth ! Upper bound of passing velocity grid normalised - real(8) :: vmint_over_vth ! Lower bound of trapped velocity grid normalised - real(8) :: vmaxt_over_vth ! Upper bound of trapped velocity grid normalised + real(dp) :: Drest(2) ! Transport coefficients (D11, D12) for trapped particles, per harmonic + real(dp) :: Tresco ! Torque from co-passing particles, per harmonic + real(dp) :: Tresctr ! Torque from counter-passing particles, per harmonic + real(dp) :: Trest ! Torque from trapped particles, per harmonic + real(dp) :: vminp_over_vth ! Lower bound of passing velocity grid normalised + real(dp) :: vmaxp_over_vth ! Upper bound of passing velocity grid normalised + real(dp) :: vmint_over_vth ! Lower bound of trapped velocity grid normalised + real(dp) :: vmaxt_over_vth ! Upper bound of trapped velocity grid normalised end type transport_harmonic_t type :: transport_data_t diff --git a/src/diag/diag_atten_map.f90 b/src/diag/diag_atten_map.f90 index e4daa944..f473d35e 100644 --- a/src/diag/diag_atten_map.f90 +++ b/src/diag/diag_atten_map.f90 @@ -1,5 +1,5 @@ module diag_atten_map - use iso_fortran_env, only: real64 + use iso_fortran_env, only: dp => real64 use fortplot, only: figure, pcolormesh, title, xlabel, ylabel, savefig use neort, only: init, check_magfie, write_magfie_data_to_files use neort_config, only: read_and_set_config @@ -22,11 +22,11 @@ subroutine run_atten_map_diag(arg_runname) character(*), intent(in) :: arg_runname logical :: file_exists integer :: i, j, nu, nm - real(real64) :: eta_max - real(real64), allocatable :: ux(:), ux_edges(:), mth_vals(:), mth_edges(:) - real(real64), allocatable :: z(:, :) - real(real64) :: v, Omth, dOmthdv, dOmthdeta, taub, bounceavg(nvar) - real(real64) :: Hmn2 + real(dp) :: eta_max + real(dp), allocatable :: ux(:), ux_edges(:), mth_vals(:), mth_edges(:) + real(dp), allocatable :: z(:, :) + real(dp) :: v, Omth, dOmthdv, dOmthdeta, taub, bounceavg(nvar) + real(dp) :: Hmn2 integer :: mth_min, mth_max integer :: u type(magfie_data_t) :: magfie_data @@ -41,30 +41,30 @@ subroutine run_atten_map_diag(arg_runname) inquire(file="plasma.in", exist=file_exists) if (file_exists) call read_and_init_plasma_input("plasma.in", s) inquire(file="profile.in", exist=file_exists) - if (file_exists) call read_and_init_profile_input("profile.in", s, R0, 1.0_real64, 1.0_real64) + if (file_exists) call read_and_init_profile_input("profile.in", s, R0, 1.0_dp, 1.0_dp) call init call check_magfie(magfie_data) call write_magfie_data_to_files(magfie_data, runname) ! Max trapped eta value consistent with earlier diagnostic - eta_max = etatp + (etadt - etatp)*(1.0_real64 - epsst_spl) + eta_max = etatp + (etadt - etatp)*(1.0_dp - epsst_spl) ! Grids nu = 80 nm = 0 - mth_min = -ceiling(2.0_real64*abs(mph*q)) - mth_max = ceiling(2.0_real64*abs(mph*q)) + mth_min = -ceiling(2.0_dp*abs(mph*q)) + mth_max = ceiling(2.0_dp*abs(mph*q)) nm = mth_max - mth_min + 1 allocate(ux(nu), ux_edges(nu+1)) allocate(mth_vals(nm), mth_edges(nm+1)) allocate(z(nm, nu)) do i = 1, nu - ux(i) = 0.2_real64 + (3.0_real64 - 0.2_real64) * real(i-1,real64) / real(nu-1,real64) + ux(i) = 0.2_dp + (3.0_dp - 0.2_dp) * real(i-1,dp) / real(nu-1,dp) end do do j = 1, nm - mth_vals(j) = real(mth_min + (j-1), real64) + mth_vals(j) = real(mth_min + (j-1), dp) end do ! Compute cell edges for pcolormesh @@ -78,9 +78,9 @@ subroutine run_atten_map_diag(arg_runname) do i = 1, nu v = ux(i)*vth call Om_th(v, eta_max, Omth, dOmthdv, dOmthdeta) - taub = 2.0_real64*acos(-1.0_real64)/abs(Omth) + taub = 2.0_dp*acos(-1.0_dp)/abs(Omth) call bounce_fast(v, eta_max, taub, bounceavg, timestep_transport) - Hmn2 = (bounceavg(3)**2 + bounceavg(4)**2)*(mi*(v*v/2.0_real64))**2 + Hmn2 = (bounceavg(3)**2 + bounceavg(4)**2)*(mi*(v*v/2.0_dp))**2 z(j, i) = nonlinear_attenuation(ux(i), eta_max, bounceavg, Omth, dOmthdv, dOmthdeta, Hmn2) end do end do @@ -110,16 +110,16 @@ subroutine run_atten_map_diag(arg_runname) end subroutine run_atten_map_diag subroutine edges_from_centers(c, e) - real(real64), intent(in) :: c(:) - real(real64), intent(out) :: e(:) + real(dp), intent(in) :: c(:) + real(dp), intent(out) :: e(:) integer :: n, i n = size(c) if (size(e) /= n+1) stop 'edges_from_centers: size mismatch' - e(1) = c(1) - 0.5_real64*(c(2)-c(1)) + e(1) = c(1) - 0.5_dp*(c(2)-c(1)) do i = 2, n - e(i) = 0.5_real64*(c(i-1)+c(i)) + e(i) = 0.5_dp*(c(i-1)+c(i)) end do - e(n+1) = c(n) + 0.5_real64*(c(n)-c(n-1)) + e(n+1) = c(n) + 0.5_dp*(c(n)-c(n-1)) end subroutine edges_from_centers end module diag_atten_map diff --git a/src/diag/diag_bounce_debug.f90 b/src/diag/diag_bounce_debug.f90 index f652be23..1f5b26c5 100644 --- a/src/diag/diag_bounce_debug.f90 +++ b/src/diag/diag_bounce_debug.f90 @@ -1,5 +1,5 @@ module diag_bounce_debug - use iso_fortran_env, only: real64 + use iso_fortran_env, only: dp => real64 use neort, only: init, check_magfie, write_magfie_data_to_files, & set_to_passing_region, set_to_trapped_region, vsteps use neort_config, only: read_and_set_config @@ -23,14 +23,14 @@ subroutine run_bounce_debug(arg_runname) character(*), intent(in) :: arg_runname logical :: file_exists integer :: i, j, mmin, mmax, nm, u - real(real64) :: vminp, vmaxp, vmint, vmaxt - real(real64) :: du, ux, v - real(real64) :: roots(100,3), eta_res(2), eta + real(dp) :: vminp, vmaxp, vmint, vmaxt + real(dp) :: du, ux, v + real(dp) :: roots(100,3), eta_res(2), eta integer :: nroots, kr - real(real64) :: Omth, dOmthdv, dOmthdeta, taub + real(dp) :: Omth, dOmthdv, dOmthdeta, taub integer :: istate integer :: istats(50) - real(real64) :: rstats(50) + real(dp) :: rstats(50) type(magfie_data_t) :: magfie_data ! Initialize environment just like main @@ -50,13 +50,13 @@ subroutine run_bounce_debug(arg_runname) call check_magfie(magfie_data) call write_magfie_data_to_files(magfie_data, runname) - vminp = 1.0d-6*vth - vmaxp = 3.0d0*vth + vminp = 1.0e-6_dp*vth + vmaxp = 3.0_dp*vth vmint = vminp vmaxt = vmaxp - mmin = -ceiling(2.0_real64*abs(mph*q)) - mmax = ceiling(2.0_real64*abs(mph*q)) + mmin = -ceiling(2.0_dp*abs(mph*q)) + mmax = ceiling(2.0_dp*abs(mph*q)) nm = mmax - mmin + 1 open(newunit=u, file=trim(arg_runname)//'_bounce_debug.txt', status='replace', action='write') @@ -68,17 +68,17 @@ subroutine run_bounce_debug(arg_runname) mth = mmin + (j - 1) ! Passing co - sign_vpar = 1.0_real64 + sign_vpar = 1.0_dp call set_to_passing_region(etamin, etamax) call scan_branch('cop', vminp, vmaxp, vsteps, u) ! Passing ctr - sign_vpar = -1.0_real64 + sign_vpar = -1.0_dp call set_to_passing_region(etamin, etamax) call scan_branch('ctr', vminp, vmaxp, vsteps, u) ! Trapped - sign_vpar = 1.0_real64 + sign_vpar = 1.0_dp call set_to_trapped_region(etamin, etamax) call scan_branch('trap', vmint, vmaxt, vsteps, u) end do @@ -87,28 +87,28 @@ end subroutine run_bounce_debug subroutine scan_branch(label, vmin, vmax, vsteps_loc, u) character(*), intent(in) :: label - real(real64), intent(in) :: vmin, vmax + real(dp), intent(in) :: vmin, vmax integer, intent(in) :: vsteps_loc integer, intent(in) :: u - real(real64) :: du, ux, v - real(real64) :: roots(100,3), eta_res(2), eta + real(dp) :: du, ux, v + real(dp) :: roots(100,3), eta_res(2), eta integer :: nroots, kr - real(real64) :: Omth, dOmthdv, dOmthdeta, taub + real(dp) :: Omth, dOmthdv, dOmthdeta, taub integer :: istate integer :: istats(50) - real(real64) :: rstats(50) + real(dp) :: rstats(50) - du = (vmax - vmin)/(real(vsteps_loc, real64)*vth) - ux = vmin/vth + du/2.0_real64 - do while (ux*vth <= vmax + 1d-12) + du = (vmax - vmin)/(real(vsteps_loc, dp)*vth) + ux = vmin/vth + du/2.0_dp + do while (ux*vth <= vmax + 1.0e-12_dp) v = ux*vth call driftorbit_coarse(v, etamin, etamax, roots, nroots) if (nroots > 0) then do kr = 1, nroots - eta_res = driftorbit_root(v, 1d-8*abs(Om_tE), roots(kr,1), roots(kr,2)) + eta_res = driftorbit_root(v, 1.0e-8_dp*abs(Om_tE), roots(kr,1), roots(kr,2)) eta = eta_res(1) call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) - taub = 2.0_real64*acos(-1.0_real64)/abs(Omth) + taub = 2.0_dp*acos(-1.0_dp)/abs(Omth) call probe_bounce(v, eta, taub, istate, istats, rstats) write(u,'(A,1X,I4,1X,F8.4,1X,ES12.5,1X,ES12.5,1X,ES12.5,1X,I4,1X,I9,1X,ES12.5)') & trim(label), mth, ux, eta, Omth, taub, istate, istats(1), rstats(6) @@ -119,26 +119,26 @@ subroutine scan_branch(label, vmin, vmax, vsteps_loc, u) end subroutine scan_branch subroutine probe_bounce(v, eta, taub, istate, istats, rstats) - real(real64), intent(in) :: v, eta, taub + real(dp), intent(in) :: v, eta, taub integer, intent(out) :: istate integer, intent(out) :: istats(:) - real(real64), intent(out) :: rstats(:) + real(dp), intent(out) :: rstats(:) - real(real64) :: t1, t2 - real(real64) :: y(nvar) - real(real64) :: atol(nvar), rtol + real(dp) :: t1, t2 + real(dp) :: y(nvar) + real(dp) :: atol(nvar), rtol integer :: neq, itask type(vode_opts) :: options - t1 = 0.0_real64 + t1 = 0.0_dp t2 = taub - y = 1.0e-15_real64 - y(1) = 0.0_real64 - y(2) = 0.0_real64 - y(3:6) = 0.0_real64 + y = 1.0e-15_dp + y(1) = 0.0_dp + y(2) = 0.0_dp + y(3:6) = 0.0_dp neq = nvar - rtol = 1.0e-8_real64 - atol = 1.0e-9_real64 + rtol = 1.0e-8_dp + atol = 1.0e-9_dp itask = 1 istate = 1 options = set_normal_opts(abserr_vector=atol, relerr=rtol) @@ -149,9 +149,9 @@ subroutine probe_bounce(v, eta, taub, istate, istats, rstats) contains subroutine timestep_wrapper(neq_, t_, y_, ydot_) integer, intent(in) :: neq_ - real(real64), intent(in) :: t_ - real(real64), intent(in) :: y_(neq_) - real(real64), intent(out) :: ydot_(neq_) + real(dp), intent(in) :: t_ + real(dp), intent(in) :: y_(neq_) + real(dp), intent(out) :: ydot_(neq_) call timestep_transport(v, eta, neq_, t_, y_, ydot_) end subroutine timestep_wrapper end subroutine probe_bounce diff --git a/src/diag/diag_bounce_nonlin.f90 b/src/diag/diag_bounce_nonlin.f90 index ea3908fd..8b8d2174 100644 --- a/src/diag/diag_bounce_nonlin.f90 +++ b/src/diag/diag_bounce_nonlin.f90 @@ -1,5 +1,5 @@ module diag_bounce_nonlin - use iso_fortran_env, only: real64 + use iso_fortran_env, only: dp => real64 use fortplot, only: figure, plot, title, xlabel, ylabel, legend, savefig, xlim, ylim use neort, only: init, check_magfie, write_magfie_data_to_files use neort_config, only: read_and_set_config @@ -22,13 +22,13 @@ module diag_bounce_nonlin subroutine run_bounce_nonlin_diag(arg_runname) character(*), intent(in) :: arg_runname logical :: file_exists - real(real64) :: v, eta_min, eta_max + real(dp) :: v, eta_min, eta_max integer :: i, j, npts - real(real64), allocatable :: eta(:) - real(real64), parameter :: ux_list(3) = [1.0_real64, 2.0_real64, 3.0_real64] - real(real64), allocatable :: att_nonlin(:, :) - real(real64) :: taub, bounceavg(nvar), Omth, dOmthdv, dOmthdeta - real(real64) :: ux, Hmn2 + real(dp), allocatable :: eta(:) + real(dp), parameter :: ux_list(3) = [1.0_dp, 2.0_dp, 3.0_dp] + real(dp), allocatable :: att_nonlin(:, :) + real(dp) :: taub, bounceavg(nvar), Omth, dOmthdv, dOmthdeta + real(dp) :: ux, Hmn2 type(magfie_data_t) :: magfie_data ! Initialize NEO-RT environment similar to neort:main @@ -42,20 +42,20 @@ subroutine run_bounce_nonlin_diag(arg_runname) if (file_exists) call read_and_init_plasma_input("plasma.in", s) inquire(file="profile.in", exist=file_exists) - if (file_exists) call read_and_init_profile_input("profile.in", s, R0, 1.0_real64, 1.0_real64) + if (file_exists) call read_and_init_profile_input("profile.in", s, R0, 1.0_dp, 1.0_dp) call init call check_magfie(magfie_data) call write_magfie_data_to_files(magfie_data, runname) ! Sample attenuation over trapped pitch range for several ux values - eta_min = etatp*(1.0_real64 + 1.0e-4_real64) - eta_max = etatp + (etadt - etatp)*(1.0_real64 - epsst_spl) + eta_min = etatp*(1.0_dp + 1.0e-4_dp) + eta_max = etatp + (etadt - etatp)*(1.0_dp - epsst_spl) npts = 100 allocate(eta(npts)) allocate(att_nonlin(npts, size(ux_list))) do i = 1, npts - eta(i) = eta_min + (eta_max-eta_min)*(real(i-1,real64)/real(npts-1,real64)) + eta(i) = eta_min + (eta_max-eta_min)*(real(i-1,dp)/real(npts-1,dp)) end do do j = 1, size(ux_list) @@ -63,9 +63,9 @@ subroutine run_bounce_nonlin_diag(arg_runname) v = ux*vth do i = 1, npts call Om_th(v, eta(i), Omth, dOmthdv, dOmthdeta) - taub = 2.0_real64*acos(-1.0_real64)/abs(Omth) + taub = 2.0_dp*acos(-1.0_dp)/abs(Omth) call bounce_fast(v, eta(i), taub, bounceavg, timestep_transport) - Hmn2 = (bounceavg(3)**2 + bounceavg(4)**2)*(mi*(v*v/2.0_real64))**2 + Hmn2 = (bounceavg(3)**2 + bounceavg(4)**2)*(mi*(v*v/2.0_dp))**2 nonlin = .true. att_nonlin(i, j) = nonlinear_attenuation(ux, eta(i), bounceavg, Omth, dOmthdv, dOmthdeta, Hmn2) @@ -87,7 +87,7 @@ subroutine run_bounce_nonlin_diag(arg_runname) call xlabel("eta") call ylabel("attenuation factor") call xlim(minval(eta), maxval(eta)) - call ylim(0.0_real64, 1.2_real64) + call ylim(0.0_dp, 1.2_dp) call legend() call savefig(trim(arg_runname)//"_bounce_nonlin.png") @@ -102,7 +102,7 @@ subroutine run_bounce_nonlin_diag(arg_runname) end subroutine run_bounce_nonlin_diag pure function to_str(x) result(s) - real(real64), intent(in) :: x + real(dp), intent(in) :: x character(len=16) :: s write(s,'(F4.1)') x s = adjustl(s) @@ -110,9 +110,9 @@ end function to_str subroutine write_data_file(fname, eta, att, uxvals) character(*), intent(in) :: fname - real(real64), intent(in) :: eta(:) - real(real64), intent(in) :: att(:, :) - real(real64), intent(in) :: uxvals(:) + real(dp), intent(in) :: eta(:) + real(dp), intent(in) :: att(:, :) + real(dp), intent(in) :: uxvals(:) integer :: i, u open(newunit=u, file=fname, status='replace', action='write') diff --git a/src/diag/diag_contrib_map.f90 b/src/diag/diag_contrib_map.f90 index 82dcb9f8..a00b93d2 100644 --- a/src/diag/diag_contrib_map.f90 +++ b/src/diag/diag_contrib_map.f90 @@ -1,5 +1,5 @@ module diag_contrib_map - use iso_fortran_env, only: real64 + use iso_fortran_env, only: dp => real64 use fortplot, only: figure, plot, pcolormesh, title, xlabel, ylabel, legend, savefig use neort, only: init, check_magfie, write_magfie_data_to_files, & set_to_passing_region, set_to_trapped_region @@ -26,18 +26,18 @@ subroutine run_contrib_diag(arg_runname) character(*), intent(in) :: arg_runname logical :: file_exists integer :: i, j, nu, nm, u - real(real64) :: v - real(real64), allocatable :: ux(:), ux_edges(:), mth_vals(:), mth_edges(:) - real(real64), allocatable :: z(:, :), y(:) + real(dp) :: v + real(dp), allocatable :: ux(:), ux_edges(:), mth_vals(:), mth_edges(:) + real(dp), allocatable :: z(:, :), y(:) integer :: mth_min, mth_max ! Temp variables for resonance calculation - real(real64) :: Omth, dOmthdv, dOmthdeta, taub, bounceavg(nvar) - real(real64) :: Hmn2, eta, du - real(real64) :: roots(100, 3) + real(dp) :: Omth, dOmthdv, dOmthdeta, taub, bounceavg(nvar) + real(dp) :: Hmn2, eta, du + real(dp) :: roots(100, 3) integer :: nroots, kr - real(real64) :: eta_res(2) - real(real64) :: contrib + real(dp) :: eta_res(2) + real(dp) :: contrib type(magfie_data_t) :: magfie_data ! Initialize environment @@ -50,7 +50,7 @@ subroutine run_contrib_diag(arg_runname) inquire(file="plasma.in", exist=file_exists) if (file_exists) call read_and_init_plasma_input("plasma.in", s) inquire(file="profile.in", exist=file_exists) - if (file_exists) call read_and_init_profile_input("profile.in", s, R0, 1.0_real64, 1.0_real64) + if (file_exists) call read_and_init_profile_input("profile.in", s, R0, 1.0_dp, 1.0_dp) call init call check_magfie(magfie_data) @@ -58,21 +58,21 @@ subroutine run_contrib_diag(arg_runname) ! Grids nu = 401 - mth_min = -ceiling(2.0_real64*abs(mph*q)) - mth_max = ceiling(2.0_real64*abs(mph*q)) + mth_min = -ceiling(2.0_dp*abs(mph*q)) + mth_max = ceiling(2.0_dp*abs(mph*q)) nm = mth_max - mth_min + 1 allocate(ux(nu), ux_edges(nu+1)) allocate(mth_vals(nm), mth_edges(nm+1)) allocate(z(nm, nu)) allocate(y(nu)) - z = 0.0_real64 - y = 0.0_real64 + z = 0.0_dp + y = 0.0_dp do i = 1, nu - ux(i) = 0.0_real64 + (4.0_real64 - 0.0_real64) * real(i-1,real64) / real(nu-1,real64) + ux(i) = 0.0_dp + (4.0_dp - 0.0_dp) * real(i-1,dp) / real(nu-1,dp) end do do j = 1, nm - mth_vals(j) = real(mth_min + (j-1), real64) + mth_vals(j) = real(mth_min + (j-1), dp) end do call edges_from_centers(ux, ux_edges) call edges_from_centers(mth_vals, mth_edges) @@ -84,18 +84,18 @@ subroutine run_contrib_diag(arg_runname) mth = nint(mth_vals(j)) do i = 1, nu v = ux(i)*vth - contrib = 0.0_real64 + contrib = 0.0_dp ! Skip exactly zero speed to avoid degenerate bounce; integrand is zero there - if (ux(i) > 0.0_real64) then + if (ux(i) > 0.0_dp) then ! Passing (co- and counter-) if enabled if (.not. nopassing) then - call accumulate_class(v, ux(i), +1.0_real64, .true., contrib) - call accumulate_class(v, ux(i), -1.0_real64, .true., contrib) + call accumulate_class(v, ux(i), +1.0_dp, .true., contrib) + call accumulate_class(v, ux(i), -1.0_dp, .true., contrib) end if ! Trapped - call accumulate_class(v, ux(i), +1.0_real64, .false., contrib) + call accumulate_class(v, ux(i), +1.0_dp, .false., contrib) end if z(j, i) = contrib @@ -142,36 +142,36 @@ subroutine run_contrib_diag(arg_runname) contains subroutine accumulate_contrib(v, ux, contrib_out) - real(real64), intent(in) :: v, ux - real(real64), intent(inout) :: contrib_out - real(real64) :: Omth, dOmthdv, dOmthdeta, taub, bounceavg(nvar), Hmn2 - real(real64) :: eta - real(real64) :: eta_res(2) - real(real64) :: roots(100, 3) + real(dp), intent(in) :: v, ux + real(dp), intent(inout) :: contrib_out + real(dp) :: Omth, dOmthdv, dOmthdeta, taub, bounceavg(nvar), Hmn2 + real(dp) :: eta + real(dp) :: eta_res(2) + real(dp) :: roots(100, 3) integer :: nroots, kr - real(real64) :: att - real(real64) :: v_eff + real(dp) :: att + real(dp) :: v_eff call driftorbit_coarse(v, etamin, etamax, roots, nroots) if (nroots == 0) return do kr = 1, nroots - eta_res = driftorbit_root(v, 1d-8*abs(Om_tE), roots(kr, 1), roots(kr, 2)) + eta_res = driftorbit_root(v, 1.0e-8_dp*abs(Om_tE), roots(kr, 1), roots(kr, 2)) eta = eta_res(1) - v_eff = max(v, 1.0e-8_real64*vth) + v_eff = max(v, 1.0e-8_dp*vth) call Om_th(v_eff, eta, Omth, dOmthdv, dOmthdeta) - taub = 2.0_real64*acos(-1.0_real64)/abs(Omth) + taub = 2.0_dp*acos(-1.0_dp)/abs(Omth) call bounce_fast(v_eff, eta, taub, bounceavg, timestep_transport) - Hmn2 = (bounceavg(3)**2 + bounceavg(4)**2)*(mi*(v_eff*v_eff/2.0_real64))**2 + Hmn2 = (bounceavg(3)**2 + bounceavg(4)**2)*(mi*(v_eff*v_eff/2.0_dp))**2 att = nonlinear_attenuation(ux, eta, bounceavg, Omth, dOmthdv, dOmthdeta, Hmn2) contrib_out = contrib_out + Tphi_int(ux, taub, Hmn2)/abs(eta_res(2)) * att end do end subroutine accumulate_contrib subroutine accumulate_class(v, ux, sign_vpar_in, is_passing, contrib_out) - real(real64), intent(in) :: v, ux, sign_vpar_in + real(dp), intent(in) :: v, ux, sign_vpar_in logical, intent(in) :: is_passing - real(real64), intent(inout) :: contrib_out + real(dp), intent(inout) :: contrib_out sign_vpar = sign_vpar_in if (is_passing) then call set_to_passing_region(etamin, etamax) @@ -184,16 +184,16 @@ end subroutine accumulate_class end subroutine run_contrib_diag subroutine edges_from_centers(c, e) - real(real64), intent(in) :: c(:) - real(real64), intent(out) :: e(:) + real(dp), intent(in) :: c(:) + real(dp), intent(out) :: e(:) integer :: n, i n = size(c) if (size(e) /= n+1) stop 'edges_from_centers: size mismatch' - e(1) = c(1) - 0.5_real64*(c(2)-c(1)) + e(1) = c(1) - 0.5_dp*(c(2)-c(1)) do i = 2, n - e(i) = 0.5_real64*(c(i-1)+c(i)) + e(i) = 0.5_dp*(c(i-1)+c(i)) end do - e(n+1) = c(n) + 0.5_real64*(c(n)-c(n-1)) + e(n+1) = c(n) + 0.5_dp*(c(n)-c(n-1)) end subroutine edges_from_centers end module diag_contrib_map diff --git a/src/do_magfie_standalone.f90 b/src/do_magfie_standalone.f90 index 25d21c62..4aab6923 100644 --- a/src/do_magfie_standalone.f90 +++ b/src/do_magfie_standalone.f90 @@ -1,34 +1,34 @@ module do_magfie_mod - + use iso_fortran_env, only: dp => real64 use util use spline implicit none - real(8), private :: s_prev = -1.0d0 - real(8), private, allocatable :: spl_val_c(:,:), spl_val_s(:,:) + real(dp), private :: s_prev = -1.0_dp + real(dp), private, allocatable :: spl_val_c(:,:), spl_val_s(:,:) ! Work arrays previously automatic on stack (size=nmode) - real(8), private, allocatable :: B0mnc(:), dB0dsmnc(:), B0mns(:), dB0dsmns(:) - real(8), private, allocatable :: costerm(:), sinterm(:) + real(dp), private, allocatable :: B0mnc(:), dB0dsmnc(:), B0mns(:), dB0dsmns(:) + real(dp), private, allocatable :: costerm(:), sinterm(:) - real(8), parameter :: sign_theta = -1.0d0 ! negative for left-handed + real(dp), parameter :: sign_theta = -1.0_dp ! negative for left-handed - real(8) :: s = 0d0, psi_pr = 0d0, Bthcov = 0d0, Bphcov = 0d0, & - dBthcovds = 0d0, dBphcovds = 0d0, & - q = 0d0, dqds = 0d0, iota = 0d0, R0 = 0d0, a = 0d0, & - eps = 0d0, B0h = 0d0, B00 = 0d0 - real(8) :: bfac = 1.0d0 + real(dp) :: s = 0.0_dp, psi_pr = 0.0_dp, Bthcov = 0.0_dp, Bphcov = 0.0_dp, & + dBthcovds = 0.0_dp, dBphcovds = 0.0_dp, & + q = 0.0_dp, dqds = 0.0_dp, iota = 0.0_dp, R0 = 0.0_dp, a = 0.0_dp, & + eps = 0.0_dp, B0h = 0.0_dp, B00 = 0.0_dp + real(dp) :: bfac = 1.0_dp ! B0h is the 0th theta harmonic of bmod on current flux surface ! and B00 the 0th theta harmonic of bmod on the innermost flux surface - real(8), allocatable, protected :: params0(:, :), modes0(:, :, :) + real(dp), allocatable, protected :: params0(:, :), modes0(:, :, :) integer, protected :: m0b, n0b, nflux, nfp, nmode - real(8), allocatable, protected :: spl_coeff1(:, :, :), spl_coeff2(:, :, :, :) + real(dp), allocatable, protected :: spl_coeff1(:, :, :), spl_coeff2(:, :, :, :) ! Work arrays for booz_to_cyl (size=nmode) - real(8), private, allocatable :: rmnc(:), rmns(:), zmnc(:), zmns(:) + real(dp), private, allocatable :: rmnc(:), rmns(:), zmnc(:), zmns(:) - real(8), parameter :: ItoB = 2.0d-1*sign_theta ! Covarient B (cgs) from I (SI) + real(dp), parameter :: ItoB = 2.0e-1_dp*sign_theta ! Covarient B (cgs) from I (SI) ! Bcov=mu0/2pi*I,mu0->4pi/c,I->10^(-1)*c*I integer :: ncol1 = 0, ncol2 = 0 ! number of columns in input file @@ -59,7 +59,7 @@ module do_magfie_mod contains subroutine set_s(s_) - real(8), intent(in) :: s_ + real(dp), intent(in) :: s_ s = s_ end subroutine set_s @@ -68,7 +68,7 @@ subroutine magfie_thread_init() ! Initialize threadprivate variables for this thread ! Must be called once per thread before using magfie routines magfie_arrays_initialized = .false. - s_prev = -1.0d0 + s_prev = -1.0_dp end subroutine magfie_thread_init subroutine read_boozer_file(path) @@ -104,15 +104,15 @@ subroutine read_boozer_file(path) end do ! Set B00 from first mode - B00 = 1.0d4*modes0(1, 1, 6)*bfac + B00 = 1.0e4_dp*modes0(1, 1, 6)*bfac end subroutine read_boozer_file subroutine init_magfie_at_s() ! Per-thread: Initialize magnetic field at current s value ! Assumes s is already set via set_s() ! Allocates threadprivate working buffers and computes s-dependent values - real(8) :: x(3), bmod, sqrtg - real(8), dimension(3) :: bder, hcovar, hctrvr, hcurl + real(dp) :: x(3), bmod, sqrtg + real(dp), dimension(3) :: bder, hcovar, hctrvr, hcurl ! Auto-initialize threadprivate state if not yet done for this thread ! (DATA/initializers don't run in worker threads) @@ -149,7 +149,7 @@ subroutine init_magfie_at_s() end if ! Initialize cache - s_prev = -1.0d0 + s_prev = -1.0_dp x(1) = s x(2) = 0.0 @@ -168,18 +168,18 @@ end subroutine do_magfie_init subroutine do_magfie(x, bmod, sqrtg, bder, hcovar, hctrvr, hcurl) ! Evaluate unperturbed axisymmetric magnetic field in point x = (s, ph, th) - real(8), dimension(:), intent(in) :: x - real(8), intent(out) :: bmod - real(8), intent(out) :: sqrtg - real(8), dimension(size(x)), intent(out) :: bder - real(8), dimension(size(x)), intent(out) :: hcovar - real(8), dimension(size(x)), intent(out) :: hctrvr - real(8), dimension(size(x)), intent(out) :: hcurl + real(dp), dimension(:), intent(in) :: x + real(dp), intent(out) :: bmod + real(dp), intent(out) :: sqrtg + real(dp), dimension(size(x)), intent(out) :: bder + real(dp), dimension(size(x)), intent(out) :: hcovar + real(dp), dimension(size(x)), intent(out) :: hctrvr + real(dp), dimension(size(x)), intent(out) :: hcurl - real(8) :: spl_val(3) - real(8) :: sqgbmod, sqgbmod2 ! sqg*B, sqg*B^2 + real(dp) :: spl_val(3) + real(dp) :: sqgbmod, sqgbmod2 ! sqg*B, sqg*B^2 - real(8) :: x1 + real(dp) :: x1 ! safety measure in order not to extrapolate ! note: this is s @@ -202,26 +202,26 @@ subroutine do_magfie(x, bmod, sqrtg, bder, hcovar, hctrvr, hcurl) ! calculate B-field from modes if (inp_swi == 8) then call cached_spline(x1, s_prev, spl_coeff2(:, :, 4, :), spl_val_c) - B0mnc(:) = 1d4*spl_val_c(1, :)*bfac - dB0dsmnc(:) = 1d4*spl_val_c(2, :)*bfac + B0mnc(:) = 1.0e4_dp*spl_val_c(1, :)*bfac + dB0dsmnc(:) = 1.0e4_dp*spl_val_c(2, :)*bfac B0h = B0mnc(1) bmod = sum(B0mnc*costerm) bder(1) = sum(dB0dsmnc*costerm)/bmod - bder(2) = 0d0 + bder(2) = 0.0_dp bder(3) = sum(-modes0(1, :, 1)*B0mnc*sinterm)/bmod else if (inp_swi == 9) then call cached_spline(x1, s_prev, spl_coeff2(:, :, 7, :), spl_val_c) - B0mnc(:) = 1d4*spl_val_c(1, :)*bfac - dB0dsmnc(:) = 1d4*spl_val_c(2, :)*bfac + B0mnc(:) = 1.0e4_dp*spl_val_c(1, :)*bfac + dB0dsmnc(:) = 1.0e4_dp*spl_val_c(2, :)*bfac call cached_spline(x1, s_prev, spl_coeff2(:, :, 8, :), spl_val_s) - B0mns(:) = 1d4*spl_val_s(1, :)*bfac - dB0dsmns(:) = 1d4*spl_val_s(2, :)*bfac + B0mns(:) = 1.0e4_dp*spl_val_s(1, :)*bfac + dB0dsmns(:) = 1.0e4_dp*spl_val_s(2, :)*bfac B0h = B0mnc(1) bmod = sum(B0mnc*costerm + B0mns*sinterm) bder(1) = sum(dB0dsmnc*costerm + dB0dsmns*sinterm)/bmod - bder(2) = 0d0 + bder(2) = 0.0_dp bder(3) = sum(-modes0(1, :, 1)*B0mnc*sinterm & + modes0(1, :, 1)*B0mns*costerm)/bmod end if @@ -230,17 +230,17 @@ subroutine do_magfie(x, bmod, sqrtg, bder, hcovar, hctrvr, hcurl) sqgbmod = sqgbmod2/bmod sqrtg = sqgbmod/bmod - hcovar(1) = 0d0 ! TODO + hcovar(1) = 0.0_dp ! TODO hcovar(2) = Bphcov/bmod hcovar(3) = Bthcov/bmod - hctrvr(1) = 0d0 + hctrvr(1) = 0.0_dp hctrvr(2) = sign_theta*psi_pr/sqgbmod hctrvr(3) = sign_theta*iota*psi_pr/sqgbmod - hcurl(1) = 0d0 ! TODO - hcurl(3) = 0d0 ! TODO - hcurl(2) = 0d0 ! TODO + hcurl(1) = 0.0_dp ! TODO + hcurl(3) = 0.0_dp ! TODO + hcurl(2) = 0.0_dp ! TODO s_prev = x1 @@ -250,7 +250,7 @@ subroutine boozer_read(filename) ! Reads Boozer in_file and converts SI to CGS integer :: ksurf, kmode - real(8) :: flux + real(dp) :: flux character(len=*) :: filename open (unit=18, file=filename, action='read', status='old') read (18, '(////)') @@ -258,7 +258,7 @@ subroutine boozer_read(filename) a = 100*a ! m -> cm R0 = 100*R0 ! m -> cm - psi_pr = 1.0d8*flux/(2*pi)*bfac ! T -> Gauss, m -> cm + psi_pr = 1.0e8_dp*flux/(2*pi)*bfac ! T -> Gauss, m -> cm nmode = (m0b + 1)*(n0b + 1) @@ -285,10 +285,10 @@ end subroutine boozer_read subroutine booz_to_cyl(x, r) - real(8), intent(in) :: x(3) ! Boozer coordinates (s, ph, th) - real(8), intent(out) :: r(3) ! Cylindrical coordinates (R, phi, Z) + real(dp), intent(in) :: x(3) ! Boozer coordinates (s, ph, th) + real(dp), intent(out) :: r(3) ! Cylindrical coordinates (R, phi, Z) - real(8) :: spl_val(3), x1 + real(dp) :: spl_val(3), x1 integer :: j @@ -299,36 +299,36 @@ subroutine booz_to_cyl(x, r) do j = 1, nmode spl_val = spline_val_0(spl_coeff2(:, :, 1, j), x1) - rmnc(j) = 1.0d2*spl_val(1) + rmnc(j) = 1.0e2_dp*spl_val(1) spl_val = spline_val_0(spl_coeff2(:, :, 2, j), x1) - rmns(j) = 1.0d2*spl_val(1) + rmns(j) = 1.0e2_dp*spl_val(1) spl_val = spline_val_0(spl_coeff2(:, :, 3, j), x1) - zmnc(j) = 1.0d2*spl_val(1) + zmnc(j) = 1.0e2_dp*spl_val(1) spl_val = spline_val_0(spl_coeff2(:, :, 4, j), x1) - zmns(j) = 1.0d2*spl_val(1) + zmns(j) = 1.0e2_dp*spl_val(1) end do r(1) = sum(rmnc*cos(modes0(1, :, 1)*x(3)) + rmns*sin(modes0(1, :, 1)*x(3))) - r(2) = 0.0d0 ! TODO: phi + r(2) = 0.0_dp ! TODO: phi r(3) = sum(zmnc*cos(modes0(1, :, 1)*x(3)) + zmns*sin(modes0(1, :, 1)*x(3))) end subroutine booz_to_cyl subroutine fast_sin_cos(m, x, sinterm_, costerm_) ! Fast sine and cosine that assumes equally spaced ascending mode numbers - real(8), intent(in) :: m(:), x - real(8), intent(out) :: sinterm_(:), costerm_(:) + real(dp), intent(in) :: m(:), x + real(dp), intent(out) :: sinterm_(:), costerm_(:) - real(8) :: dm - complex(8) :: fourier_factor, rotation + real(dp) :: dm + complex(dp) :: fourier_factor, rotation integer :: j dm = m(2) - m(1) fourier_factor = exp(imun*m(1)*x) rotation = exp(imun*dm*x) - costerm_ = (0.0d0, 0.0d0) - sinterm_ = (0.0d0, 0.0d0) + costerm_ = (0.0_dp, 0.0_dp) + sinterm_ = (0.0_dp, 0.0_dp) do j = 1, size(m) costerm_(j) = real(fourier_factor) sinterm_(j) = imag(fourier_factor) @@ -339,27 +339,27 @@ end subroutine fast_sin_cos end module do_magfie_mod module do_magfie_pert_mod - + use iso_fortran_env, only: dp => real64 use util use spline use do_magfie_mod, only: s, bfac, inp_swi implicit none - real(8), private :: s_prev = -1.0d0 - real(8), private, allocatable :: spl_val_c(:,:), spl_val_s(:,:) + real(dp), private :: s_prev = -1.0_dp + real(dp), private, allocatable :: spl_val_c(:,:), spl_val_s(:,:) - real(8), allocatable, protected :: params(:, :), modes(:, :, :) + real(dp), allocatable, protected :: params(:, :), modes(:, :, :) integer, protected :: mb, nb, nflux, nfp, nmode - real(8), allocatable, protected :: spl_coeff1(:, :, :), spl_coeff2(:, :, :, :) + real(dp), allocatable, protected :: spl_coeff1(:, :, :), spl_coeff2(:, :, :, :) ! Work arrays (size=nmode) - real(8), private, allocatable :: Bmnc(:), Bmns(:) + real(dp), private, allocatable :: Bmnc(:), Bmns(:) integer :: ncol1, ncol2 ! number of columns in input file - real(8) :: mph ! toroidal perturbation mode (threadprivate) - real(8) :: mph_shared = 0d0 ! shared copy for namelist input (when pertfile=.false.) + real(dp) :: mph ! toroidal perturbation mode (threadprivate) + real(dp) :: mph_shared = 0.0_dp ! shared copy for namelist input (when pertfile=.false.) ! Initialization flag for threadprivate allocatable arrays logical, save :: magfie_pert_arrays_initialized = .false. @@ -384,7 +384,7 @@ subroutine magfie_pert_thread_init() ! Initialize threadprivate variables for this thread ! Must be called once per thread before using magfie_pert routines magfie_pert_arrays_initialized = .false. - s_prev = -1.0d0 + s_prev = -1.0_dp end subroutine magfie_pert_thread_init subroutine read_boozer_pert_file(path) @@ -427,8 +427,8 @@ subroutine init_magfie_pert_at_s() ! Per-thread: Initialize perturbation field at current s value ! Assumes s is already set via set_s() ! Allocates threadprivate working buffers and computes s-dependent values - real(8) :: x(3) - complex(8) :: dummy + real(dp) :: x(3) + complex(dp) :: dummy ! Auto-initialize threadprivate state if not yet done for this thread if (magfie_pert_thread_init_state /= MAGFIE_PERT_INIT_SENTINEL) then @@ -453,7 +453,7 @@ subroutine init_magfie_pert_at_s() end if ! Initialize cache - s_prev = -1.0d0 + s_prev = -1.0_dp ! Compute mph at current s mph = nfp*modes(1, 1, 2) @@ -465,7 +465,7 @@ subroutine init_magfie_pert_at_s() end subroutine init_magfie_pert_at_s subroutine set_mph(mph_value) - real(8), intent(in) :: mph_value + real(dp), intent(in) :: mph_value mph = mph_value mph_shared = mph_value end subroutine set_mph @@ -484,10 +484,10 @@ subroutine do_magfie_pert_init(path) end subroutine do_magfie_pert_init subroutine do_magfie_pert_amp(x, bamp) - real(8), dimension(:), intent(in) :: x - complex(8), intent(out) :: bamp + real(dp), dimension(:), intent(in) :: x + complex(dp), intent(out) :: bamp - real(8) :: x1 + real(dp) :: x1 ! safety measure in order not to extrapolate x1 = max(params(1, 1), x(1)) @@ -496,13 +496,13 @@ subroutine do_magfie_pert_amp(x, bamp) ! calculate B-field from modes if (inp_swi == 8) then call cached_spline(x1, s_prev, spl_coeff2(:, :, 4, :), spl_val_c) - Bmnc(:) = 1d4*spl_val_c(1, :)*bfac + Bmnc(:) = 1.0e4_dp*spl_val_c(1, :)*bfac bamp = sum(Bmnc*cos(modes(1, :, 1)*x(3))) else if (inp_swi == 9) then call cached_spline(x1, s_prev, spl_coeff2(:, :, 7, :), spl_val_c) call cached_spline(x1, s_prev, spl_coeff2(:, :, 8, :), spl_val_s) - Bmnc(:) = 1d4*spl_val_c(1, :)*bfac - Bmns(:) = 1d4*spl_val_s(1, :)*bfac + Bmnc(:) = 1.0e4_dp*spl_val_c(1, :)*bfac + Bmns(:) = 1.0e4_dp*spl_val_s(1, :)*bfac bamp = fast_fourier_sum(Bmnc, Bmns, modes(1, :, 1), x(3)) end if @@ -510,9 +510,9 @@ subroutine do_magfie_pert_amp(x, bamp) end subroutine do_magfie_pert_amp subroutine do_magfie_pert(x, bmod) - real(8), dimension(:), intent(in) :: x - real(8), intent(out) :: bmod - complex(8) :: bamp + real(dp), dimension(:), intent(in) :: x + real(dp), intent(out) :: bmod + complex(dp) :: bamp call do_magfie_pert_amp(x, bamp) bmod = real(sum(bamp*exp(imun*nfp*modes(1, :, 2)*x(2)))) @@ -521,7 +521,7 @@ end subroutine do_magfie_pert subroutine boozer_read_pert(filename) integer :: ksurf, kmode - real(8) :: flux, dummy + real(dp) :: flux, dummy character(len=*) :: filename open (unit=18, file=filename) read (18, '(////)') @@ -550,11 +550,11 @@ subroutine boozer_read_pert(filename) end subroutine boozer_read_pert subroutine check_equal_space_ascending(m, kmode) - real(8), dimension(:), intent(in) :: m + real(dp), dimension(:), intent(in) :: m integer, intent(in) :: kmode - real(8), parameter :: tol = 1.0d-10 - real(8) :: diff, diff2 + real(dp), parameter :: tol = 1.0e-10_dp + real(dp) :: diff, diff2 if (kmode < 3) return @@ -566,10 +566,10 @@ subroutine check_equal_space_ascending(m, kmode) end subroutine check_equal_space_ascending subroutine check_equal(m, kmode) - real(8), dimension(:), intent(in) :: m + real(dp), dimension(:), intent(in) :: m integer, intent(in) :: kmode - real(8), parameter :: tol = 1.0d-10 + real(dp), parameter :: tol = 1.0e-10_dp if (kmode < 2) return @@ -580,19 +580,19 @@ end subroutine check_equal function fast_fourier_sum(fmnc, fmns, m, x) ! Fast Fourier sum that assumes equally spaced ascending mode numbers - real(8), dimension(:), intent(in) :: fmnc, fmns, m - real(8), intent(in) :: x - complex(8) :: fast_fourier_sum + real(dp), dimension(:), intent(in) :: fmnc, fmns, m + real(dp), intent(in) :: x + complex(dp) :: fast_fourier_sum - real(8) :: dm - complex(8) :: fourier_factor, rotation + real(dp) :: dm + complex(dp) :: fourier_factor, rotation integer :: j dm = m(2) - m(1) fourier_factor = exp(imun*m(1)*x) rotation = exp(imun*dm*x) - fast_fourier_sum = (0.0d0, 0.0d0) + fast_fourier_sum = (0.0_dp, 0.0_dp) do j = 1, size(m) fast_fourier_sum = fast_fourier_sum + & (fmnc(j) - imun*fmns(j))*fourier_factor diff --git a/src/driftorbit.f90 b/src/driftorbit.f90 index 584a0445..3e7c47d0 100644 --- a/src/driftorbit.f90 +++ b/src/driftorbit.f90 @@ -3,6 +3,7 @@ ! Christopher Albert, since 2015 module driftorbit + use iso_fortran_env, only: dp => real64 use util use do_magfie_mod use do_magfie_pert_mod, only: do_magfie_pert_amp, mph @@ -13,11 +14,11 @@ module driftorbit implicit none ! Normalization factor for radial electric field - real(8) :: efac = 1d0 + real(dp) :: efac = 1.0_dp ! Harmonics TODO: make dynamic, multiple harmonics ! Default values are overridden by config file in driftorbit_test:read_control - real(8) :: epsmn = 1d0 ! perturbation amplitude B1/B0 + real(dp) :: epsmn = 1.0_dp ! perturbation amplitude B1/B0 integer :: m0 = 1 ! Boozer poloidal perturbation mode integer :: mth = 1 ! canonical poloidal mode logical :: magdrift = .true. ! consider magnetic drift @@ -26,22 +27,22 @@ module driftorbit logical :: comptorque = .true. ! compute torque ! Flux surface TODO: make a dynamic, multiple flux surfaces support - real(8) :: dVds = 0d0, etadt = 0d0, etatp = 0d0 - real(8) :: etamin = 0d0, etamax = 0d0 + real(dp) :: dVds = 0.0_dp, etadt = 0.0_dp, etatp = 0.0_dp + real(dp) :: etamin = 0.0_dp, etamax = 0.0_dp ! TODO: better B0 calculation (average magnetic field on flux surface) - real(8) :: B0 = 0d0 - real(8) :: Bmin = 0d0, Bmax = 0d0 + real(dp) :: B0 = 0.0_dp + real(dp) :: Bmin = 0.0_dp, Bmax = 0.0_dp - real(8), parameter :: epst_spl = 1d-6, epsp_spl = 1d-6 ! dist to tpb for spline - real(8), parameter :: epsst_spl = 1d-3, epssp_spl = 1d-3 ! dist to deep for spline - real(8), parameter :: epst = 1d-8, epsp = 1d-8 ! smallest eta distance to tp bound + real(dp), parameter :: epst_spl = 1.0e-6_dp, epsp_spl = 1.0e-6_dp ! dist to tpb for spline + real(dp), parameter :: epsst_spl = 1.0e-3_dp, epssp_spl = 1.0e-3_dp ! dist to deep for spline + real(dp), parameter :: epst = 1.0e-8_dp, epsp = 1.0e-8_dp ! smallest eta distance to tp bound ! Number of levels for coarse root finding integer, parameter :: nlev = 100 - real(8) :: sign_vpar = 1d0, sign_vpar_htheta = 1d0 + real(dp) :: sign_vpar = 1.0_dp, sign_vpar_htheta = 1.0_dp ! Nonlinear calculation switch logical :: nonlin = .false. diff --git a/src/freq.f90 b/src/freq.f90 index 95b6ecf2..e4a9c0f0 100644 --- a/src/freq.f90 +++ b/src/freq.f90 @@ -1,4 +1,5 @@ module neort_freq + use iso_fortran_env, only: dp => real64 use logger, only: debug, trace, get_log_level, LOG_TRACE use util, only: pi use spline, only: spline_coeff, spline_val_0 @@ -12,18 +13,18 @@ module neort_freq ! For splining in the trapped eta region integer, parameter :: netaspl = 100 - real(8), allocatable :: OmtB_spl_coeff(:, :) - real(8), allocatable :: Omth_spl_coeff(:, :) - real(8), allocatable :: vres_spl_coeff(:, :) + real(dp), allocatable :: OmtB_spl_coeff(:, :) + real(dp), allocatable :: Omth_spl_coeff(:, :) + real(dp), allocatable :: vres_spl_coeff(:, :) ! For splining in the passing eta region integer, parameter :: netaspl_pass = 100 - real(8), allocatable :: OmtB_pass_spl_coeff(:, :) - real(8), allocatable :: Omth_pass_spl_coeff(:, :) - real(8), allocatable :: vres_pass_spl_coeff(:, :) + real(dp), allocatable :: OmtB_pass_spl_coeff(:, :) + real(dp), allocatable :: Omth_pass_spl_coeff(:, :) + real(dp), allocatable :: vres_pass_spl_coeff(:, :) - real(8) :: k_taub_p=0d0, d_taub_p=0d0, k_taub_t=0d0, d_taub_t=0d0 ! extrapolation at tp bound - real(8) :: k_OmtB_p=0d0, d_Omtb_p=0d0, k_Omtb_t=0d0, d_Omtb_t=0d0 ! extrapolation at tp bound + real(dp) :: k_taub_p=0.0_dp, d_taub_p=0.0_dp, k_taub_t=0.0_dp, d_taub_t=0.0_dp ! extrapolation at tp bound + real(dp) :: k_OmtB_p=0.0_dp, d_Omtb_p=0.0_dp, k_Omtb_t=0.0_dp, d_Omtb_t=0.0_dp ! extrapolation at tp bound ! Initialization flags for threadprivate allocatable arrays logical :: freq_trapped_initialized = .false. @@ -49,24 +50,24 @@ subroutine freq_thread_init() call vode_thread_init() freq_trapped_initialized = .false. freq_passing_initialized = .false. - k_taub_p = 0d0 - d_taub_p = 0d0 - k_taub_t = 0d0 - d_taub_t = 0d0 - k_OmtB_p = 0d0 - d_Omtb_p = 0d0 - k_Omtb_t = 0d0 - d_Omtb_t = 0d0 + k_taub_p = 0.0_dp + d_taub_p = 0.0_dp + k_taub_t = 0.0_dp + d_taub_t = 0.0_dp + k_OmtB_p = 0.0_dp + d_Omtb_p = 0.0_dp + k_Omtb_t = 0.0_dp + d_Omtb_t = 0.0_dp end subroutine freq_thread_init subroutine init_canon_freq_trapped_spline ! Initialise splines for canonical frequencies of trapped orbits - real(8) :: etarange(netaspl), Om_tB_v(netaspl), Omth_v(netaspl) + real(dp) :: etarange(netaspl), Om_tB_v(netaspl), Omth_v(netaspl) integer :: k - real(8) :: aa, b - real(8) :: taub0, taub1, leta0, leta1, OmtB0, OmtB1 - real(8) :: v, eta, taub, taub_est, bounceavg(nvar) + real(dp) :: aa, b + real(dp) :: taub0, taub1, leta0, leta1, OmtB0, OmtB1 + real(dp) :: v, eta, taub, taub_est, bounceavg(nvar) ! Auto-initialize threadprivate state if not yet done for this thread if (freq_thread_init_state /= FREQ_INIT_SENTINEL) then @@ -76,16 +77,16 @@ subroutine init_canon_freq_trapped_spline call trace('init_canon_freq_trapped_spline') - taub0 = 0d0 - taub1 = 0d0 - leta0 = 0d0 - leta1 = 0d0 - OmtB0 = 0d0 - OmtB1 = 0d0 + taub0 = 0.0_dp + taub1 = 0.0_dp + leta0 = 0.0_dp + leta1 = 0.0_dp + OmtB0 = 0.0_dp + OmtB1 = 0.0_dp v = vth - etamin = (1d0 + epst)*etatp - etamax = etatp + (etadt - etatp)*(1d0 - epsst_spl) + etamin = (1.0_dp + epst)*etatp + etamax = etatp + (etadt - etatp)*(1.0_dp - epsst_spl) ! Allocate coefficient arrays for trapped region splines (safe for undefined allocation status) if (.not. freq_trapped_initialized) then if (allocated(Omth_spl_coeff)) deallocate(Omth_spl_coeff) @@ -108,10 +109,10 @@ subroutine init_canon_freq_trapped_spline ! logspace for eta b = log(epst_spl) - aa = 1d0/(netaspl - 1d0)*(log(etamax/etamin - 1d0) - b) + aa = 1.0_dp/(netaspl - 1.0_dp)*(log(etamax/etamin - 1.0_dp) - b) do k = netaspl - 1, 0, -1 - eta = etamin*(1d0 + exp(aa*k + b)) + eta = etamin*(1.0_dp + exp(aa*k + b)) etarange(k + 1) = eta if (get_log_level() >= LOG_TRACE) then write(*,'(A,I4,A,ES12.5)') '[TRACE] init_canon_freq_trapped_spline k=', k, ' eta=', eta @@ -157,11 +158,11 @@ end subroutine init_canon_freq_trapped_spline subroutine init_canon_freq_passing_spline ! Initialise splines for canonical frequencies of passing orbits - real(8) :: etarange(netaspl_pass), Om_tB_v(netaspl_pass), Omth_v(netaspl_pass) - real(8) :: aa, b + real(dp) :: etarange(netaspl_pass), Om_tB_v(netaspl_pass), Omth_v(netaspl_pass) + real(dp) :: aa, b integer :: k - real(8) :: leta0, leta1, taub0, taub1, OmtB0, OmtB1 - real(8) :: v, eta, taub, taub_est, bounceavg(nvar) + real(dp) :: leta0, leta1, taub0, taub1, OmtB0, OmtB1 + real(dp) :: v, eta, taub, taub_est, bounceavg(nvar) ! Auto-initialize threadprivate state if not yet done for this thread if (freq_thread_init_state /= FREQ_INIT_SENTINEL) then @@ -171,12 +172,12 @@ subroutine init_canon_freq_passing_spline call trace('init_canon_freq_passing_spline') - taub0 = 0d0 - taub1 = 0d0 - leta0 = 0d0 - leta1 = 0d0 - OmtB0 = 0d0 - OmtB1 = 0d0 + taub0 = 0.0_dp + taub1 = 0.0_dp + leta0 = 0.0_dp + leta1 = 0.0_dp + OmtB0 = 0.0_dp + OmtB1 = 0.0_dp v = vth etamin = etatp*epssp_spl @@ -202,10 +203,10 @@ subroutine init_canon_freq_passing_spline end if b = log((etamax - etamin)/etamax) - aa = 1d0/(netaspl_pass - 1d0)*(log(epsp_spl) - b) + aa = 1.0_dp/(netaspl_pass - 1.0_dp)*(log(epsp_spl) - b) do k = netaspl_pass - 1, 0, -1 - eta = etamax*(1d0 - exp(aa*k + b)) + eta = etamax*(1.0_dp - exp(aa*k + b)) etarange(k + 1) = eta if (get_log_level() >= LOG_TRACE) then write(*,'(A,I4,A,ES12.5)') '[TRACE] init_canon_freq_passing_spline k=', k, ' eta=', eta @@ -249,10 +250,10 @@ end subroutine init_canon_freq_passing_spline subroutine Om_tB(v, eta, OmtB, dOmtBdv, dOmtBdeta) ! returns bounce averaged toroidal magnetic drift frequency ! and derivatives w.r.t. v and eta - real(8), intent(in) :: v, eta - real(8), intent(out) :: OmtB, dOmtBdv, dOmtBdeta - real(8) :: splineval(3) - real(8) :: Omth, dOmthdv, dOmthdeta + real(dp), intent(in) :: v, eta + real(dp), intent(out) :: OmtB, dOmtBdv, dOmtBdeta + real(dp) :: splineval(3) + real(dp) :: Omth, dOmthdv, dOmthdeta if (eta > etatp) then if (eta > etatp*(1 + epst_spl)) then splineval = spline_val_0(OmtB_spl_coeff, eta) @@ -273,22 +274,22 @@ subroutine Om_tB(v, eta, OmtB, dOmtBdv, dOmtBdeta) end if end if OmtB = splineval(1)*v**2 - dOmtBdv = 2d0*splineval(1)*v + dOmtBdv = 2.0_dp*splineval(1)*v dOmtBdeta = splineval(2)*v**2 end subroutine Om_tB subroutine Om_ph(v, eta, Omph, dOmphdv, dOmphdeta) ! returns canonical toroidal frequency ! and derivatives w.r.t. v and eta - real(8), intent(in) :: v, eta - real(8), intent(out) :: Omph, dOmphdv, dOmphdeta - real(8) :: Omth, dOmthdv, dOmthdeta - real(8) :: OmtB, dOmtBdv, dOmtBdeta + real(dp), intent(in) :: v, eta + real(dp), intent(out) :: Omph, dOmphdv, dOmphdeta + real(dp) :: Omth, dOmthdv, dOmthdeta + real(dp) :: OmtB, dOmtBdv, dOmtBdeta if (eta > etatp) then Omph = Om_tE - dOmphdv = 0d0 - dOmphdeta = 0d0 + dOmphdv = 0.0_dp + dOmphdeta = 0.0_dp if (magdrift) then call Om_tB(v, eta, OmtB, dOmtBdv, dOmtBdeta) Omph = Omph + OmtB @@ -312,23 +313,23 @@ end subroutine Om_ph subroutine Om_th(v, eta, Omth, dOmthdv, dOmthdeta) ! returns canonical poloidal frequency ! and derivatives w.r.t. v and eta - real(8), intent(in) :: v, eta - real(8), intent(out) :: Omth, dOmthdv, dOmthdeta - real(8) :: splineval(3) + real(dp), intent(in) :: v, eta + real(dp), intent(out) :: Omth, dOmthdv, dOmthdeta + real(dp) :: splineval(3) if (eta > etatp) then if (eta > etatp*(1 + epst_spl)) then splineval = spline_val_0(Omth_spl_coeff, eta) else ! extrapolation - splineval(1) = 2d0*pi/(k_taub_t*log(eta - etatp) + d_taub_t) - splineval(2) = -splineval(1)**2/(2d0*pi)*k_taub_t/(eta - etatp) + splineval(1) = 2.0_dp*pi/(k_taub_t*log(eta - etatp) + d_taub_t) + splineval(2) = -splineval(1)**2/(2.0_dp*pi)*k_taub_t/(eta - etatp) end if else if (eta < etatp*(1 - epsp_spl)) then splineval = spline_val_0(Omth_pass_spl_coeff, eta) else ! extrapolation - splineval(1) = 2d0*pi/(k_taub_p*log(etatp - eta) + d_taub_p) - splineval(2) = -splineval(1)**2/(2d0*pi)*k_taub_p/(eta - etatp) + splineval(1) = 2.0_dp*pi/(k_taub_p*log(etatp - eta) + d_taub_p) + splineval(2) = -splineval(1)**2/(2.0_dp*pi)*k_taub_p/(eta - etatp) end if end if Omth = sign_vpar*splineval(1)*v @@ -337,22 +338,22 @@ subroutine Om_th(v, eta, Omth, dOmthdv, dOmthdeta) end subroutine Om_th subroutine d_Om_ds(v, eta, taub_estimate, dOmthds, dOmphds) - real(8), intent(in) :: v, eta, taub_estimate - real(8), intent(out) :: dOmthds, dOmphds - real(8) :: s0, ds, bounceavg(nvar) - real(8) :: taub, taub_est, Omth, Omph_noE + real(dp), intent(in) :: v, eta, taub_estimate + real(dp), intent(out) :: dOmthds, dOmphds + real(dp) :: s0, ds, bounceavg(nvar) + real(dp) :: taub, taub_est, Omth, Omph_noE call trace('d_Om_ds') ! store current flux surface values s0 = s - ds = 2d-8 - s = s0 - ds/2d0 + ds = 2.0e-8_dp + s = s0 - ds/2.0_dp taub_est = bounce_time(v, eta, taub_estimate) taub = taub_est call bounce_fast(v, eta, taub, bounceavg, timestep) - Omth = sign_vpar_htheta*2d0*pi/taub + Omth = sign_vpar_htheta*2.0_dp*pi/taub if (magdrift) then if (eta > etatp) then Omph_noE = bounceavg(3)*v**2 @@ -361,27 +362,27 @@ subroutine d_Om_ds(v, eta, taub_estimate, dOmthds, dOmphds) end if else if (eta > etatp) then - Omph_noE = 0d0 + Omph_noE = 0.0_dp else Omph_noE = Omth/iota end if end if - s = s0 + ds/2d0 + s = s0 + ds/2.0_dp taub_est = bounce_time(v, eta, taub_estimate) taub = taub_est call bounce_fast(v, eta, taub, bounceavg, timestep) - dOmthds = sign_vpar_htheta*(2d0*pi/taub - sign_vpar_htheta*Omth)/ds + dOmthds = sign_vpar_htheta*(2.0_dp*pi/taub - sign_vpar_htheta*Omth)/ds if (magdrift) then if (eta > etatp) then dOmphds = dOm_tEds + (bounceavg(3)*v**2 - Omph_noE)/ds else - dOmphds = dOm_tEds + (bounceavg(3)*v**2 + (2d0*pi/taub)/iota - Omph_noE)/ds + dOmphds = dOm_tEds + (bounceavg(3)*v**2 + (2.0_dp*pi/taub)/iota - Omph_noE)/ds end if else if (eta > etatp) then dOmphds = dOm_tEds else - dOmphds = dOm_tEds + ((2d0*pi/taub)/iota - Omph_noE)/ds + dOmphds = dOm_tEds + ((2.0_dp*pi/taub)/iota - Omph_noE)/ds end if end if diff --git a/src/magfie.f90 b/src/magfie.f90 index 41b2a5ad..87f92c9a 100644 --- a/src/magfie.f90 +++ b/src/magfie.f90 @@ -1,4 +1,5 @@ module neort_magfie + use iso_fortran_env, only: dp => real64 use util, only: disp, pi use do_magfie_mod, only: do_magfie, eps, iota use neort_orbit, only: th0 @@ -12,11 +13,11 @@ module neort_magfie subroutine init_flux_surface_average(s) ! Calculate the flux surface areas for normalization ! Note: thread initialization is now handled by auto-init pattern in each module - real(8), intent(in) :: s + real(dp), intent(in) :: s integer, parameter :: nth = 1000 integer :: k - real(8) :: thrange(nth), dth - real(8) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) + real(dp) :: thrange(nth), dth + real(dp) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) character(len=256) :: buffer write(buffer, "(A,ES12.5)") " s: ", s @@ -26,17 +27,17 @@ subroutine init_flux_surface_average(s) dth = thrange(2) - thrange(1) x(1) = s - x(2) = 0d0 - x(3) = 0d0 + x(2) = 0.0_dp + x(3) = 0.0_dp - dVds = 0d0 - B0 = 0d0 + dVds = 0.0_dp + B0 = 0.0_dp write(buffer, "(A,ES12.5)") " eps orig: ", eps call log_result(buffer) - eps = 0d0 + eps = 0.0_dp - Bmin = -1d0 - Bmax = 0d0 + Bmin = -1.0_dp + Bmax = 0.0_dp do k = 1, nth x(3) = thrange(k) @@ -53,12 +54,12 @@ subroutine init_flux_surface_average(s) if (bmod > Bmax) Bmax = bmod end do - dVds = 2d0*pi*dVds - B0 = B0/(2d0*pi) + dVds = 2.0_dp*pi*dVds + B0 = B0/(2.0_dp*pi) eps = eps/(B0*pi) - etatp = 1d0/Bmax - etadt = 1d0/Bmin + etatp = 1.0_dp/Bmax + etadt = 1.0_dp/Bmin write(buffer, "(A,ES12.5)") " eps calc: ", eps call log_result(buffer) diff --git a/src/neort.f90 b/src/neort.f90 index 44770002..3531e065 100644 --- a/src/neort.f90 +++ b/src/neort.f90 @@ -1,4 +1,5 @@ module neort + use iso_fortran_env, only: dp => real64 use logger, only: debug, set_log_level, get_log_level, log_result, LOG_INFO use neort_datatypes, only: magfie_data_t, transport_data_t, transport_harmonic_t use neort_profiles, only: init_thermodynamic_forces, init_profiles, vth, dvthds, ni1, & @@ -32,7 +33,7 @@ end subroutine init pure subroutine set_to_trapped_region(eta_min, eta_max) use driftorbit, only: epst, etatp, etadt - real(8), intent(out) :: eta_min, eta_max + real(dp), intent(out) :: eta_min, eta_max eta_min = (1 + epst)*etatp eta_max = (1 - epst)*etadt @@ -41,35 +42,37 @@ end subroutine set_to_trapped_region pure subroutine set_to_passing_region(eta_min, eta_max) use driftorbit, only: epsp, etatp - real(8), intent(out) :: eta_min, eta_max + real(dp), intent(out) :: eta_min, eta_max eta_min = epsp*etatp eta_max = (1 - epsp)*etatp end subroutine set_to_passing_region subroutine check_magfie(data) - use do_magfie_mod, only: do_magfie, sign_theta - use do_magfie_pert_mod, only: do_magfie_pert_amp - use driftorbit + use util, only: pi, c, imun, eV + use do_magfie_mod, only: do_magfie, sign_theta, R0, a, psi_pr, iota, q, Bthcov, Bphcov, & + dBthcovds, dBphcovds, eps + use do_magfie_pert_mod, only: do_magfie_pert_amp, mph + use driftorbit, only: B0, etatp, etadt, M_t, Om_tE, m0, pertfile, nonlin, epsmn, dVds type(magfie_data_t), intent(out) :: data integer, parameter :: nth = 50 integer :: k - real(8) :: thmin, thmax - real(8) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) - real(8) :: Dp, Drp - real(8) :: dpp, dhh, fpeff - complex(8) :: bn + real(dp) :: thmin, thmax + real(dp) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) + real(dp) :: D_plateau, Drp + real(dp) :: dpp, dhh, fpeff + complex(dp) :: bn - Dp = pi*vth**3/(16d0*R0*iota*(qi*B0/(mi*c))**2) + D_plateau = pi*vth**3/(16.0_dp*R0*iota*(qi*B0/(mi*c))**2) Drp = 4*mph*q/(eps**2*sqrt(pi)) ! actually this is Drp/Dp thmin = -pi thmax = pi x(1) = s - x(2) = 0d0 - x(3) = 0d0 + x(2) = 0.0_dp + x(3) = 0.0_dp data%params%s = s data%params%R0 = R0 @@ -86,20 +89,20 @@ subroutine check_magfie(data) data%params%dVds = dVds data%params%M_t = M_t data%params%Om_tE = Om_tE - data%params%Om_tBref = c*mi*vth**2/(2d0*qi*sign_theta*psi_pr) + data%params%Om_tBref = c*mi*vth**2/(2.0_dp*qi*sign_theta*psi_pr) data%params%vth = vth - data%params%T_in_eV = mi/2d0*vth**2/eV - data%params%m0 = 1d0*m0 - data%params%n0 = 1d0*mph - data%params%Dp = Dp + data%params%T_in_eV = mi/2.0_dp*vth**2/eV + data%params%m0 = 1.0_dp*m0 + data%params%n0 = 1.0_dp*mph + data%params%Dp = D_plateau data%params%Drp_over_Dp = Drp data%params%etatp = etatp data%params%etadt = etadt data%params%pertfile = pertfile data%params%nonlin = nonlin - dpp = 0d0 - dhh = 0d0 - fpeff = 0d0 + dpp = 0.0_dp + dhh = 0.0_dp + fpeff = 0.0_dp data%params%dpp = dpp data%params%dhh = dhh data%params%fpeff = fpeff @@ -130,29 +133,32 @@ subroutine check_magfie(data) end subroutine check_magfie subroutine compute_transport(result_) - use driftorbit + use driftorbit, only: mth, M_t, R0, etamin, etamax, sign_vpar, nopassing, comptorque, & + dVds, mph, dOm_tEds, dM_tds + use neort_profiles, only: Om_tE + use do_magfie_mod, only: q type(transport_data_t), intent(out) :: result_ integer :: j, idx - ! Transport coefficients D11, D12 in approximate effective radius r=2d0/a*sqrt(s) - real(8) :: Dco(2), Dctr(2), Dt(2) + ! Transport coefficients D11, D12 in approximate effective radius r=2.0_dp/a*sqrt(s) + real(dp) :: Dco(2), Dctr(2), Dt(2) ! Torque density dTphi_int/ds for integration over normalized toroidal flux s - real(8) :: Tco, Tctr, Tt + real(dp) :: Tco, Tctr, Tt integer :: mthmin, mthmax call debug('compute_transport') if (allocated(result_%harmonics)) deallocate(result_%harmonics) - Dco = 0d0 - Dctr = 0d0 - Dt = 0d0 - Tco = 0d0 - Tctr = 0d0 - Tt = 0d0 + Dco = 0.0_dp + Dctr = 0.0_dp + Dt = 0.0_dp + Tco = 0.0_dp + Tctr = 0.0_dp + Tt = 0.0_dp Om_tE = vth*M_t/R0 dOm_tEds = vth*dM_tds/R0 + M_t*dvthds/R0 @@ -188,25 +194,25 @@ subroutine compute_transport(result_) result_%torque%Tctr = Tctr result_%torque%Tt = Tt else - result_%torque%s = 0d0 - result_%torque%dVds = 0d0 - result_%torque%M_t = 0d0 - result_%torque%Tco = 0d0 - result_%torque%Tctr = 0d0 - result_%torque%Tt = 0d0 + result_%torque%s = 0.0_dp + result_%torque%dVds = 0.0_dp + result_%torque%M_t = 0.0_dp + result_%torque%Tco = 0.0_dp + result_%torque%Tctr = 0.0_dp + result_%torque%Tt = 0.0_dp end if call debug('compute_transport complete') end subroutine compute_transport subroutine compute_transport_harmonic(j, Dco, Dctr, Dt, Tco, Tctr, Tt, harmonic) - use driftorbit + use driftorbit, only: mth, M_t, etamin, etamax, sign_vpar, nopassing integer, intent(in) :: j - real(8), intent(inout) :: Dco(2), Dctr(2), Dt(2), Tco, Tctr, Tt + real(dp), intent(inout) :: Dco(2), Dctr(2), Dt(2), Tco, Tctr, Tt type(transport_harmonic_t), intent(out) :: harmonic - real(8) :: Dresco(2), Dresctr(2), Drest(2), Tresco, Tresctr, Trest - real(8) :: vminp, vmaxp, vmint, vmaxt + real(dp) :: Dresco(2), Dresctr(2), Drest(2), Tresco, Tresctr, Trest + real(dp) :: vminp, vmaxp, vmint, vmaxt character(len=256) :: buffer write(buffer, '(A,ES12.5,A,I0)') "compute_transport_harmonic: M_t = ", M_t, ", mth = ", j @@ -214,15 +220,15 @@ subroutine compute_transport_harmonic(j, Dco, Dctr, Dt, Tco, Tctr, Tt, harmonic) mth = j - Drest = 0d0 - Dresco = 0d0 - Dresctr = 0d0 - Tresco = 0d0 - Tresctr = 0d0 - Trest = 0d0 + Drest = 0.0_dp + Dresco = 0.0_dp + Dresctr = 0.0_dp + Tresco = 0.0_dp + Tresctr = 0.0_dp + Trest = 0.0_dp - vminp = 1d-6*vth - vmaxp = 3d0*vth + vminp = 1.0e-6_dp*vth + vmaxp = 3.0_dp*vth vmint = vminp vmaxt = vmaxp @@ -297,7 +303,7 @@ subroutine write_magfie_data_to_files(data, base_path) write (unit, *) "check_magfie: R0 = ", data%params%R0 write (unit, *) "check_magfie: a = ", data%params%a write (unit, *) "check_magfie: eps = ", data%params%eps - write (unit, *) "check_magfie: A = ", 1.0d0 / data%params%eps + write (unit, *) "check_magfie: A = ", 1.0_dp / data%params%eps write (unit, *) "check_magfie: psi_pr = ", data%params%psi_pr write (unit, *) "check_magfie: B0 = ", data%params%B0 write (unit, *) "check_magfie: Bthcov = ", data%params%Bthcov @@ -349,7 +355,7 @@ subroutine write_transport_data_to_files(data, base_path) character(len=*), intent(in) :: base_path integer :: k - real(8) :: total_D1, total_D2 + real(dp) :: total_D1, total_D2 integer, parameter :: unit1 = 9 integer, parameter :: unit2 = 10 diff --git a/src/nonlin.f90 b/src/nonlin.f90 index b729957e..7a20bf0d 100644 --- a/src/nonlin.f90 +++ b/src/nonlin.f90 @@ -1,4 +1,5 @@ module neort_nonlin + use iso_fortran_env, only: dp => real64 use logger, only: debug, trace, get_log_level, LOG_TRACE use util, only: c, mi, pi use collis_alp, only: coleff @@ -12,16 +13,16 @@ module neort_nonlin contains function nonlinear_attenuation(ux, eta, bounceavg, Omth, dOmthdv, dOmthdeta, Hmn2) - real(8), intent(in) :: ux, eta, bounceavg(nvar), Omth, dOmthdv, dOmthdeta, Hmn2 - real(8) :: nonlinear_attenuation + real(dp), intent(in) :: ux, eta, bounceavg(nvar), Omth, dOmthdv, dOmthdeta, Hmn2 + real(dp) :: nonlinear_attenuation - real(8) :: dpp, dhh, fpeff, dres, dnorm, Omph, dOmphdv, dOmphdeta, dOmdv, & + real(dp) :: dpp, dhh, fpeff, dres, dnorm, Omph, dOmphdv, dOmphdeta, dOmdv, & dOmdeta, Ompr, dOmphds, dOmthds, dOmdpph, v call trace('nonlinear_attenuation') if (.not. nonlin) then - nonlinear_attenuation = 1d0 + nonlinear_attenuation = 1.0_dp return end if @@ -29,7 +30,7 @@ function nonlinear_attenuation(ux, eta, bounceavg, Omth, dOmthdv, dOmthdeta, Hmn if (nonlin) then call Om_ph(v, eta, Omph, dOmphdv, dOmphdeta) - call d_Om_ds(v, eta, 2d0*pi/Omth, dOmthds, dOmphds) + call d_Om_ds(v, eta, 2.0_dp*pi/Omth, dOmthds, dOmphds) dOmdv = mth*dOmthdv + mph*dOmphdv dOmdeta = mth*dOmthdeta + mph*dOmphdeta dOmdpph = -(qi/c*iota*sign_theta*psi_pr)**(-1)*(mth*dOmthds + mph*dOmphds) @@ -38,7 +39,7 @@ function nonlinear_attenuation(ux, eta, bounceavg, Omth, dOmthdv, dOmthdeta, Hmn dhh = vth*dhh dpp = vth**3*dpp dres = dpp*(dOmdv/Ompr)**2 + dhh*eta*(bounceavg(5) - eta)*(dOmdeta/Ompr)**2 - dnorm = dres*sqrt(abs(Ompr))/sqrt(abs(Hmn2))**(3d0/2d0) + dnorm = dres*sqrt(abs(Ompr))/sqrt(abs(Hmn2))**(3.0_dp/2.0_dp) call attenuation_factor(dnorm, nonlinear_attenuation) end if @@ -46,17 +47,17 @@ function nonlinear_attenuation(ux, eta, bounceavg, Omth, dOmthdv, dOmthdeta, Hmn end function nonlinear_attenuation pure function omega_prime(ux, eta, bounceavg, Omth, dOmdv, dOmdeta, dOmdpph) - real(8), intent(in) :: ux, eta, bounceavg(nvar), Omth, dOmdv, dOmdeta, dOmdpph - real(8) :: omega_prime + real(dp), intent(in) :: ux, eta, bounceavg(nvar), Omth, dOmdv, dOmdeta, dOmdpph + real(dp) :: omega_prime - real(8) :: ma, mb, mc, md, me, mf, dvdJ, detadJ + real(dp) :: ma, mb, mc, md, me, mf, dvdJ, detadJ ma = mi*(ux*vth)*mi*c/qi*eta mb = mi*(ux*vth)**2/2*mi*c/qi - mc = mi/(2d0*Omth)*(ux*vth)*(1d0 - eta*bounceavg(6)) - md = mi*(ux*vth)**2/2d0*Omth + mc = mi/(2.0_dp*Omth)*(ux*vth)*(1.0_dp - eta*bounceavg(6)) + md = mi*(ux*vth)**2/2.0_dp*Omth me = -mth/mph - mf = 1d0/mph + mf = 1.0_dp/mph dvdJ = mb*me/(ma*md*mf - mb*mc*mf) detadJ = ma*me/(mb*mc*mf - ma*md*mf) diff --git a/src/orbit.f90 b/src/orbit.f90 index 1b5ef292..9f2cf66b 100644 --- a/src/orbit.f90 +++ b/src/orbit.f90 @@ -1,4 +1,5 @@ module neort_orbit + use iso_fortran_env, only: dp => real64 use logger, only: debug, trace, get_log_level, LOG_TRACE, error use util, only: imun, pi, mi, qi, c use spline, only: spline_coeff, spline_val_0 @@ -13,7 +14,7 @@ module neort_orbit implicit none integer, parameter :: nvar = 7 - real(8) :: th0 = 0d0 + real(dp) :: th0 = 0.0_dp logical :: noshear = .false. ! neglect magnetic shear @@ -21,11 +22,12 @@ module neort_orbit interface subroutine timestep_i(v, eta, neq, t, y, ydot) - real(8), intent(in) :: v, eta + import :: dp + real(dp), intent(in) :: v, eta integer, intent(in) :: neq - real(8), intent(in) :: t - real(8), intent(in) :: y(neq) - real(8), intent(out) :: ydot(neq) + real(dp), intent(in) :: t + real(dp), intent(in) :: y(neq) + real(dp), intent(out) :: ydot(neq) end subroutine timestep_i end interface @@ -36,7 +38,7 @@ subroutine dvode_error_context(where, v_in, eta_in, tcur, tout, ist) use driftorbit, only: etatp, etadt, etamin, etamax, mth, mph, sign_vpar use neort_profiles, only: vth, Om_tE character(*), intent(in) :: where - real(8), intent(in) :: v_in, eta_in, tcur, tout + real(dp), intent(in) :: v_in, eta_in, tcur, tout integer, intent(in) :: ist character(len=512) :: msg character(len=64) :: reg @@ -59,7 +61,7 @@ subroutine dvode_error_context(where, v_in, eta_in, tcur, tout, ist) end subroutine dvode_error_context pure function to_es(x) result(sout) - real(8), intent(in) :: x + real(dp), intent(in) :: x character(len=24) :: sout write(sout,'(ES12.5)') x end function to_es @@ -73,48 +75,48 @@ end function to_i subroutine bounce(v, eta, taub, bounceavg, taub_estimate) ! calculate all bounce averages - real(8), intent(in) :: v, eta - real(8), intent(out) :: taub, bounceavg(nvar) - real(8), optional :: taub_estimate ! estimated bounce time (user input) - real(8) :: findroot_res(nvar + 1) - real(8) :: bmod, htheta - real(8) :: y0(nvar) + real(dp), intent(in) :: v, eta + real(dp), intent(out) :: taub, bounceavg(nvar) + real(dp), optional :: taub_estimate ! estimated bounce time (user input) + real(dp) :: findroot_res(nvar + 1) + real(dp) :: bmod, htheta + real(dp) :: y0(nvar) ! Initialize bounce-averated quantities y0. Their meaning ! is defined inside subroutine timestep (thin orbit integration) call evaluate_bfield_local(bmod, htheta) - sign_vpar_htheta = sign(1d0, htheta)*sign_vpar - y0 = 1d-15 + sign_vpar_htheta = sign(1.0_dp, htheta)*sign_vpar + y0 = 1.0e-15_dp y0(1) = th0 ! poloidal angle theta y0(2) = sign_vpar_htheta*vpar(v, eta, bmod) ! parallel velocity vpar - y0(3) = 0d0 ! toroidal velocity v_ph for drift frequency Om_ph - y0(4) = 0d0 ! perturbed Hamiltonian real part - y0(5) = 0d0 ! perturbed Hamiltonian imaginary part - y0(6) = 0d0 ! 1/abs(B) - ! y0(7) = 0d0 ! abs(B) + y0(3) = 0.0_dp ! toroidal velocity v_ph for drift frequency Om_ph + y0(4) = 0.0_dp ! perturbed Hamiltonian real part + y0(5) = 0.0_dp ! perturbed Hamiltonian imaginary part + y0(6) = 0.0_dp ! 1/abs(B) + ! y0(7) = 0.0_dp ! abs(B) ! If bounce time estimate exists (elliptic integrals), ! initialize taub with it, owtherwise estimate here. if (present(taub_estimate)) then taub = taub_estimate else - taub = 2.0*pi/abs(vperp(v, eta, bmod)*iota/R0*sqrt(eps/2d0)) + taub = 2.0*pi/abs(vperp(v, eta, bmod)*iota/R0*sqrt(eps/2.0_dp)) end if ! Look for exactly one orbit turn via root-finding. ! Start by looking for 5 points per turn. - findroot_res = bounce_integral(v, eta, nvar, y0, taub/5d0, timestep) + findroot_res = bounce_integral(v, eta, nvar, y0, taub/5.0_dp, timestep) taub = findroot_res(1) bounceavg = findroot_res(2:)/taub end subroutine bounce subroutine evaluate_bfield_local(bmod, htheta) - real(8), intent(out) :: bmod, htheta - real(8) :: sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) + real(dp), intent(out) :: bmod, htheta + real(dp) :: sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) x(1) = s - x(2) = 0d0 + x(2) = 0.0_dp x(3) = th0 call do_magfie(x, bmod, sqrtg, hder, hcovar, hctrvr, hcurl) htheta = hctrvr(3) @@ -122,53 +124,53 @@ end subroutine evaluate_bfield_local pure function vpar(v, eta, bmod) ! parallel velocity - real(8) :: vpar - real(8), intent(in) :: v, eta, bmod - vpar = v*sqrt(1d0 - eta*bmod) + real(dp) :: vpar + real(dp), intent(in) :: v, eta, bmod + vpar = v*sqrt(1.0_dp - eta*bmod) if (isnan(vpar)) then - vpar = 0d0 + vpar = 0.0_dp end if end function vpar pure function vperp(v, eta, bmod) ! perpendicular velocity - real(8) :: vperp - real(8), intent(in) :: v, eta, bmod + real(dp) :: vperp + real(dp), intent(in) :: v, eta, bmod vperp = v*sqrt(eta*bmod) if (isnan(vperp)) then - vperp = 0d0 + vperp = 0.0_dp end if end function vperp subroutine bounce_fast(v, eta, taub, bounceavg, ts, istate_out) use dvode_f90_m - real(8), intent(in) :: v, eta, taub - real(8), intent(out) :: bounceavg(nvar) + real(dp), intent(in) :: v, eta, taub + real(dp), intent(out) :: bounceavg(nvar) procedure(timestep_i) :: ts integer, intent(out), optional :: istate_out - real(8) :: t1, t2, bmod, htheta - real(8) :: y(nvar) - real(8) :: atol(nvar), rtol + real(dp) :: t1, t2, bmod, htheta + real(dp) :: y(nvar) + real(dp) :: atol(nvar), rtol integer :: neq, itask, istate type(vode_opts) :: options call trace('bounce_fast') - t1 = 0d0 + t1 = 0.0_dp t2 = taub call evaluate_bfield_local(bmod, htheta) - sign_vpar_htheta = sign(1d0, htheta)*sign_vpar - y = 1d-15 + sign_vpar_htheta = sign(1.0_dp, htheta)*sign_vpar + y = 1.0e-15_dp y(1) = th0 y(2) = sign_vpar_htheta*vpar(v, eta, bmod) - y(3:6) = 0d0 + y(3:6) = 0.0_dp neq = nvar - rtol = 1d-9 - atol = 1d-10 + rtol = 1.0e-9_dp + atol = 1.0e-10_dp itask = 1 istate = 1 @@ -188,9 +190,9 @@ subroutine bounce_fast(v, eta, taub, bounceavg, ts, istate_out) subroutine timestep_wrapper(neq_, t_, y_, ydot_) ! Wrapper routine for timestep to work with VODE integer, intent(in) :: neq_ - real(8), intent(in) :: t_ - real(8), intent(in) :: y_(neq_) - real(8), intent(out) :: ydot_(neq_) + real(dp), intent(in) :: t_ + real(dp), intent(in) :: y_(neq_) + real(dp), intent(out) :: ydot_(neq_) call ts(v, eta, neq_, t_, y_, ydot_) end subroutine timestep_wrapper @@ -198,17 +200,17 @@ end subroutine bounce_fast function bounce_time(v, eta, taub_estimate) result(taub) - real(8), intent(in) :: v, eta - real(8), intent(in), optional :: taub_estimate - real(8) :: taub + real(dp), intent(in) :: v, eta + real(dp), intent(in), optional :: taub_estimate + real(dp) :: taub integer, parameter :: neq = 2 - real(8) :: y0(neq), roots(neq+1) - real(8) :: bmod, htheta + real(dp) :: y0(neq), roots(neq+1) + real(dp) :: bmod, htheta call trace('bounce_time') call evaluate_bfield_local(bmod, htheta) - sign_vpar_htheta = sign(1d0, htheta)*sign_vpar + sign_vpar_htheta = sign(1.0_dp, htheta)*sign_vpar y0(1) = th0 ! poloidal angle theta y0(2) = sign_vpar_htheta*vpar(v, eta, bmod) ! parallel velocity vpar @@ -216,7 +218,7 @@ function bounce_time(v, eta, taub_estimate) result(taub) if (present(taub_estimate)) then taub = taub_estimate else - taub = 2.0*pi/abs(vperp(v, eta, bmod)*iota/R0*sqrt(eps/2d0)) + taub = 2.0*pi/abs(vperp(v, eta, bmod)*iota/R0*sqrt(eps/2.0_dp)) end if roots = bounce_integral(v, eta, neq, y0, taub, timestep_poloidal_motion) @@ -225,15 +227,15 @@ function bounce_time(v, eta, taub_estimate) result(taub) end function bounce_time subroutine timestep_poloidal_motion(v, eta, neq, t, y, ydot) - real(8), intent(in) :: v, eta + real(dp), intent(in) :: v, eta integer, intent(in) :: neq - real(8), intent(in) :: t, y(neq) - real(8), intent(out) :: ydot(neq) + real(dp), intent(in) :: t, y(neq) + real(dp), intent(out) :: ydot(neq) - real(8) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) + real(dp) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) x(1) = s - x(2) = 0d0 + x(2) = 0.0_dp x(3) = y(1) call do_magfie(x, bmod, sqrtg, hder, hcovar, hctrvr, hcurl) @@ -241,12 +243,12 @@ subroutine timestep_poloidal_motion(v, eta, neq, t, y, ydot) end subroutine timestep_poloidal_motion pure subroutine poloidal_velocity(v, eta, bmod, hthctr, hderth, v_par, ydot) - real(8), intent(in) :: v, eta, bmod, hthctr, hderth - real(8), intent(in) :: v_par - real(8), intent(out) :: ydot(2) + real(dp), intent(in) :: v, eta, bmod, hthctr, hderth + real(dp), intent(in) :: v_par + real(dp), intent(out) :: ydot(2) ydot(1) = v_par*hthctr ! theta - ydot(2) = -v**2*eta/2d0*hthctr*hderth*bmod ! v_par + ydot(2) = -v**2*eta/2.0_dp*hthctr*hderth*bmod ! v_par end subroutine poloidal_velocity function bounce_integral(v, eta, neq, y0, dt, ts) @@ -255,26 +257,26 @@ function bounce_integral(v, eta, neq, y0, dt, ts) ! use dvode_f90_m - real(8) :: bounce_integral(neq + 1) - real(8), intent(in) :: v, eta + real(dp) :: bounce_integral(neq + 1) + real(dp), intent(in) :: v, eta integer, intent(in) :: neq - real(8), intent(in) :: y0(neq), dt + real(dp), intent(in) :: y0(neq), dt procedure(timestep_i) :: ts integer :: n integer :: k, state, rootstate - real(8) :: ti, told - real(8) :: y(neq), yold(neq) + real(dp) :: ti, told + real(dp) :: y(neq), yold(neq) logical :: passing - real(8) :: atol(neq), rtol, tout + real(dp) :: atol(neq), rtol, tout integer :: itask, istate type(vode_opts) :: options - rtol = 1d-9 - atol = 1d-10 + rtol = 1.0e-9_dp + atol = 1.0e-10_dp itask = 1 istate = 1 @@ -289,7 +291,7 @@ function bounce_integral(v, eta, neq, y0, dt, ts) y = y0 yold = y0 - ti = 0d0 + ti = 0.0_dp state = 1 if (get_log_level() >= LOG_TRACE) then write(*,'(A,2ES12.5,2A)') '[TRACE] bounce_integral start v,eta=', v, eta, ' pass=', merge('T','F',eta T) end associate GOUT(1) = sign_vpar_htheta*(Y(1) - th0) ! trapped orbit return to starting point - GOUT(2) = sign_vpar_htheta*(2d0*pi - (Y(1) - th0)) ! passing orbit return + GOUT(2) = sign_vpar_htheta*(2.0_dp*pi - (Y(1) - th0)) ! passing orbit return return end subroutine bounceroots @@ -363,19 +365,19 @@ subroutine timestep(v, eta, neq, t, y, ydot) ! More integrands may be added starting from y(3) ! - real(8), intent(in) :: v, eta + real(dp), intent(in) :: v, eta integer, intent(in) :: neq - real(8), intent(in) :: t - real(8), intent(in) :: y(neq) - real(8), intent(out) :: ydot(neq) + real(dp), intent(in) :: t + real(dp), intent(in) :: y(neq) + real(dp), intent(out) :: ydot(neq) - real(8) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) - real(8) :: Om_tB_v - real(8) :: shearterm + real(dp) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) + real(dp) :: Om_tB_v + real(dp) :: shearterm x(1) = s - x(2) = 0d0 + x(2) = 0.0_dp x(3) = y(1) call do_magfie(x, bmod, sqrtg, hder, hcovar, hctrvr, hcurl) @@ -384,15 +386,15 @@ subroutine timestep(v, eta, neq, t, y, ydot) shearterm = 0 end if - Om_tB_v = mi*c*q/(2d0*qi*sign_theta*psi_pr*bmod)*( & ! Om_tB/v**2 - -(2d0 - eta*bmod)*bmod*hder(1) & - + 2d0*(1d0 - eta*bmod)*hctrvr(3)* & + Om_tB_v = mi*c*q/(2.0_dp*qi*sign_theta*psi_pr*bmod)*( & ! Om_tB/v**2 + -(2.0_dp - eta*bmod)*bmod*hder(1) & + + 2.0_dp*(1.0_dp - eta*bmod)*hctrvr(3)* & (dBthcovds + q*dBphcovds + shearterm)) ydot(1) = y(2)*hctrvr(3) ! theta - ydot(2) = -0.5d0*v**2*eta*hctrvr(3)*hder(3)*bmod ! v_par + ydot(2) = -0.5_dp*v**2*eta*hctrvr(3)*hder(3)*bmod ! v_par ydot(3) = Om_tB_v ! for bounce average of Om_tB/v**2 - ydot(4:) = 0d0 ! remaining integrands not computed here + ydot(4:) = 0.0_dp ! remaining integrands not computed here end subroutine timestep end module neort_orbit diff --git a/src/profiles.f90 b/src/profiles.f90 index 9a64d979..fd046cfc 100644 --- a/src/profiles.f90 +++ b/src/profiles.f90 @@ -8,21 +8,21 @@ module neort_profiles implicit none ! Thermal velocity, Mach number, and their derivatives over s_tor - real(dp) :: vth = 0d0, dvthds = 0d0, M_t = 0d0, dM_tds = 0d0 + real(dp) :: vth = 0.0_dp, dvthds = 0.0_dp, M_t = 0.0_dp, dM_tds = 0.0_dp ! Electric precession frequency and its derivative over s_tor - real(8) :: Om_tE = 0d0, dOm_tEds = 0d0 + real(dp) :: Om_tE = 0.0_dp, dOm_tEds = 0.0_dp ! density and temperature profiles and their derivatives over s_tor - real(dp) :: ni1 = 0d0, ni2 = 0d0, Ti1 = 0d0, Ti2 = 0d0, Te = 0d0, dni1ds = 0d0, & - dni2ds = 0d0, dTi1ds = 0d0, dTi2ds = 0d0, dTeds = 0d0 + real(dp) :: ni1 = 0.0_dp, ni2 = 0.0_dp, Ti1 = 0.0_dp, Ti2 = 0.0_dp, Te = 0.0_dp, dni1ds = 0.0_dp, & + dni2ds = 0.0_dp, dTi1ds = 0.0_dp, dTi2ds = 0.0_dp, dTeds = 0.0_dp ! Thermodynamic forces in radial variable s_tor - real(dp) :: A1 = 0d0, A2 = 0d0 + real(dp) :: A1 = 0.0_dp, A2 = 0.0_dp ! Plasma profile spline data integer :: nplasma_global = 0 - real(dp) :: am1_global = 0d0, am2_global = 0d0, Z1_global = 0d0, Z2_global = 0d0 + real(dp) :: am1_global = 0.0_dp, am2_global = 0.0_dp, Z1_global = 0.0_dp, Z2_global = 0.0_dp real(dp), allocatable :: plasma_spl_coeff(:,:,:) ! Rotation profile spline data @@ -38,7 +38,7 @@ subroutine init_profiles(R0) real(dp), intent(in) :: R0 Om_tE = vth * M_t / R0 ! toroidal ExB drift frequency - dOm_tEds = 0d0 + dOm_tEds = 0.0_dp end subroutine init_profiles subroutine prepare_plasma_splines(nplasma, am1, am2, Z1, Z2, plasma) @@ -71,7 +71,7 @@ subroutine init_plasma_at_s() ! Assumes s is already set via set_s() and prepare_plasma_splines() was called use do_magfie_mod, only: s real(dp) :: spl_val(3) - real(dp), parameter :: pmass = 1.6726d-24 + real(dp), parameter :: pmass = 1.6726e-24_dp real(dp) :: amb, Zb, dchichi, slowrate, dchichi_norm, slowrate_norm real(dp) :: v0, ebeam @@ -98,14 +98,14 @@ subroutine init_plasma_at_s() ! Compute derived quantities qi = Z1_global*qe mi = am1_global*mu - vth = sqrt(2d0*Ti1*ev/mi) - dvthds = 0.5d0*sqrt(2d0*ev/(mi*Ti1))*dTi1ds + vth = sqrt(2.0_dp*Ti1*ev/mi) + dvthds = 0.5_dp*sqrt(2.0_dp*ev/(mi*Ti1))*dTi1ds ! Call collision routine v0 = vth - amb = 2d0 - Zb = 1d0 - ebeam = amb*pmass*v0**2/(2d0*ev) + amb = 2.0_dp + Zb = 1.0_dp + ebeam = amb*pmass*v0**2/(2.0_dp*ev) call loacol_nbi(amb, am1_global, am2_global, Zb, Z1_global, Z2_global, & ni1, ni2, Ti1, Ti2, Te, ebeam, v0, & dchichi, slowrate, dchichi_norm, slowrate_norm) @@ -158,7 +158,7 @@ end subroutine read_and_init_plasma_input subroutine prepare_profile_splines(data) ! Main thread only: Prepare shared spline coefficients from rotation profile data ! This should be called ONCE before parallel region - real(8), intent(in) :: data(:, :) + real(dp), intent(in) :: data(:, :) ! Deallocate if already allocated if (allocated(Mt_spl_coeff)) deallocate(Mt_spl_coeff) @@ -172,8 +172,8 @@ subroutine init_profile_at_s(R0, efac, bfac) ! Per-thread: Interpolate rotation profile at current s value ! Assumes s is already set and prepare_profile_splines() was called use do_magfie_mod, only: s - real(8), intent(in) :: R0, efac, bfac - real(8) :: splineval(3) + real(dp), intent(in) :: R0, efac, bfac + real(dp) :: splineval(3) ! Interpolate M_t at s (writes threadprivate variables) splineval = spline_val_0(Mt_spl_coeff, s) @@ -191,9 +191,9 @@ subroutine read_and_init_profile_input(path, s_in, R0, efac, bfac) ! This combines readdata, prepare_profile_splines, and init_profile_at_s use do_magfie_mod, only: s character(len=*), intent(in) :: path - real(8), intent(in) :: s_in, R0, efac, bfac + real(dp), intent(in) :: s_in, R0, efac, bfac - real(8), allocatable :: data(:, :) + real(dp), allocatable :: data(:, :) call readdata(path, 2, data) ! allocates data @@ -210,7 +210,7 @@ subroutine init_thermodynamic_forces(psi_pr, q) real(dp), intent(in) :: psi_pr ! toroidal flux at plasma boundary == dpsi_tor/ds real(dp), intent(in) :: q ! safety factor - A1 = dni1ds/ni1 - qi/(Ti1*ev)*sign_theta*psi_pr/(q*c)*Om_tE - 3d0/2d0*dTi1ds/Ti1 + A1 = dni1ds/ni1 - qi/(Ti1*ev)*sign_theta*psi_pr/(q*c)*Om_tE - 3.0_dp/2.0_dp*dTi1ds/Ti1 A2 = dTi1ds/Ti1 end subroutine init_thermodynamic_forces diff --git a/src/resonance.f90 b/src/resonance.f90 index 3a25577f..4603bf39 100644 --- a/src/resonance.f90 +++ b/src/resonance.f90 @@ -1,24 +1,25 @@ module neort_resonance + use iso_fortran_env, only: dp => real64 use neort_freq, only: Om_th, Om_ph, d_Om_ds use driftorbit, only: mth, mph, nlev, vth, sign_vpar implicit none contains subroutine driftorbit_coarse(v, eta_min, eta_max, roots, nroots) - real(8), intent(in) :: v, eta_min, eta_max - real(8), intent(out) :: roots(:, :) + real(dp), intent(in) :: v, eta_min, eta_max + real(dp), intent(out) :: roots(:, :) integer, intent(out) :: nroots - real(8) :: deta - real(8) :: Omph, dOmphdv, dOmphdeta - real(8) :: Omth, dOmthdv, dOmthdeta - real(8) :: res, dresdv, dresdeta - real(8) :: resold, dresdvold, dresdetaold - real(8) :: eta + real(dp) :: deta + real(dp) :: Omph, dOmphdv, dOmphdeta + real(dp) :: Omth, dOmthdv, dOmthdeta + real(dp) :: res, dresdv, dresdeta + real(dp) :: resold, dresdvold, dresdetaold + real(dp) :: eta integer :: k, ninterv ninterv = size(roots, 1) - deta = (eta_max - eta_min)*1d0/ninterv + deta = (eta_max - eta_min)*1.0_dp/ninterv nroots = 0 do k = 0, ninterv @@ -29,7 +30,7 @@ subroutine driftorbit_coarse(v, eta_min, eta_max, roots, nroots) dresdv = mth*dOmthdv + mph*dOmphdv dresdeta = mth*dOmthdeta + mph*dOmphdeta if (k > 0) then - if (sign(1d0, res) /= sign(1d0, resold)) then + if (sign(1.0_dp, res) /= sign(1.0_dp, resold)) then nroots = nroots + 1 roots(nroots, 1) = eta - deta roots(nroots, 2) = eta @@ -43,9 +44,9 @@ end subroutine driftorbit_coarse function driftorbit_nroot(v, eta_min, eta_max) integer :: driftorbit_nroot - real(8), intent(in) :: v - real(8), intent(in) :: eta_min, eta_max - real(8) :: roots(nlev, 3) + real(dp), intent(in) :: v + real(dp), intent(in) :: eta_min, eta_max + real(dp) :: roots(nlev, 3) call driftorbit_coarse(v, eta_min, eta_max, roots, driftorbit_nroot) end function driftorbit_nroot @@ -53,16 +54,16 @@ end function driftorbit_nroot function driftorbit_root(v, tol, eta_min, eta_max) use logger, only: warning - real(8) :: driftorbit_root(2) - real(8), intent(in) :: v, tol, eta_min, eta_max - real(8) :: res, res_old, eta0, eta_old - real(8) :: Omph, dOmphdv, dOmphdeta - real(8) :: Omth, dOmthdv, dOmthdeta + real(dp) :: driftorbit_root(2) + real(dp), intent(in) :: v, tol, eta_min, eta_max + real(dp) :: res, res_old, eta0, eta_old + real(dp) :: Omph, dOmphdv, dOmphdeta + real(dp) :: Omth, dOmthdv, dOmthdeta integer :: maxit, k, state - real(8) :: etamin2, etamax2 + real(dp) :: etamin2, etamax2 logical :: slope_pos - real(8) :: resmin, resmax - real(8) :: eta + real(dp) :: resmin, resmax + real(dp) :: eta character(len=1024) :: msg character, parameter :: TAB = char(9) @@ -71,8 +72,8 @@ function driftorbit_root(v, tol, eta_min, eta_max) maxit = 100 state = -2 eta0 = eta - eta_old = 0d0 - res = 0d0 + eta_old = 0.0_dp + res = 0.0_dp etamin2 = eta_min etamax2 = eta_max @@ -117,11 +118,11 @@ function driftorbit_root(v, tol, eta_min, eta_max) ((.not. slope_pos) .and. res < 0)) then etamax2 = eta eta_old = eta - eta = (eta + etamin2)/2d0 + eta = (eta + etamin2)/2.0_dp else etamin2 = eta eta_old = eta - eta = (eta + etamax2)/2d0 + eta = (eta + etamax2)/2.0_dp end if end do if (state < 0) then diff --git a/src/transport.f90 b/src/transport.f90 index d652d91b..92599bb3 100644 --- a/src/transport.f90 +++ b/src/transport.f90 @@ -1,4 +1,5 @@ module neort_transport + use iso_fortran_env, only: dp => real64 use util, only: imun, pi, c, qi use logger, only: trace, debug, warning, error use do_magfie_mod, only: do_magfie, s, a, R0, iota, q, psi_pr, eps, & @@ -16,7 +17,7 @@ module neort_transport implicit none - real(8) :: Omth, dOmthdv, dOmthdeta + real(dp) :: Omth, dOmthdv, dOmthdeta contains @@ -24,8 +25,8 @@ pure function fmt_dbg(msg1, v1, msg2, v2, msg3, v3, msg4, v4) result(s) ! Helper to compose a short debug line character(*), intent(in) :: msg1, msg2 character(*), intent(in), optional :: msg3, msg4 - real(8), intent(in) :: v1, v2 - real(8), intent(in), optional :: v3, v4 + real(dp), intent(in) :: v1, v2 + real(dp), intent(in), optional :: v3, v4 character(len=256) :: s character(len=64) :: a1, a2, a3, a4 a3 = ''; a4 = '' @@ -45,71 +46,71 @@ end function fmt_dbg ! original contains follows pure function D11int(ux, taub, Hmn2) - real(8) :: D11int - real(8), intent(in) :: ux, taub, Hmn2 + real(dp) :: D11int + real(dp), intent(in) :: ux, taub, Hmn2 - D11int = pi**(3d0/2d0)*mph**2*c**2*q*vth & + D11int = pi**(3.0_dp/2.0_dp)*mph**2*c**2*q*vth & /(qi**2*dVds*abs(psi_pr))*ux**3*exp(-ux**2)*taub*Hmn2 end function D11int pure function D12int(ux, taub, Hmn2) - real(8) :: D12int - real(8), intent(in) :: ux, taub, Hmn2 + real(dp) :: D12int + real(dp), intent(in) :: ux, taub, Hmn2 D12int = D11int(ux, taub, Hmn2)*ux**2 end function D12int pure function Tphi_int(ux, taub, Hmn2) - real(8) :: Tphi_int - real(8), intent(in) :: ux, taub, Hmn2 + real(dp) :: Tphi_int + real(dp), intent(in) :: ux, taub, Hmn2 - Tphi_int = sign(1d0, psi_pr*q*sign_theta)*pi**(3d0/2d0)*mph**2*ni1*c*vth/qi & + Tphi_int = sign(1.0_dp, psi_pr*q*sign_theta)*pi**(3.0_dp/2.0_dp)*mph**2*ni1*c*vth/qi & *ux**3*exp(-ux**2)*taub*Hmn2*(A1 + A2*ux**2) end function Tphi_int subroutine compute_transport_integral(vmin, vmax, vsteps, D, T) ! compute transport integral via midpoint rule - real(8), intent(in) :: vmin, vmax + real(dp), intent(in) :: vmin, vmax integer, intent(in) :: vsteps - real(8), intent(out) :: D(2), T ! Transport coefficients D and torque density T - real(8) :: Dp, dsdreff ! Plateau diffusion coefficient and ds/dreff=<|grad s|> - real(8) :: ux, du, dD11, dD12, dT, v, eta - real(8) :: eta_res(2) - real(8) :: taub, bounceavg(nvar) + real(dp), intent(out) :: D(2), T ! Transport coefficients D and torque density T + real(dp) :: D_plateau, dsdreff ! Plateau diffusion coefficient and ds/dreff=<|grad s|> + real(dp) :: ux, du, dD11, dD12, dT, v, eta + real(dp) :: eta_res(2) + real(dp) :: taub, bounceavg(nvar) integer :: istate_dv - real(8) :: Hmn2, attenuation_factor - real(8) :: roots(nlev, 3) + real(dp) :: Hmn2, attenuation_factor + real(dp) :: roots(nlev, 3) integer :: nroots, kr, ku call debug(fmt_dbg('compute_transport_integral: vmin=', vmin, ' vmax=', vmax, ' vsteps=', dble(vsteps))) - D = 0d0 - T = 0d0 + D = 0.0_dp + T = 0.0_dp du = (vmax - vmin)/(vsteps*vth) - ux = vmin/vth + du/2d0 + ux = vmin/vth + du/2.0_dp do ku = 1, vsteps v = ux*vth call driftorbit_coarse(v, etamin, etamax, roots, nroots) if (nroots == 0) continue do kr = 1, nroots - eta_res = driftorbit_root(v, 1d-8*abs(Om_tE), roots(kr, 1), roots(kr, 2)) + eta_res = driftorbit_root(v, 1.0e-8_dp*abs(Om_tE), roots(kr, 1), roots(kr, 2)) eta = eta_res(1) call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) - taub = 2d0*pi/abs(Omth) + taub = 2.0_dp*pi/abs(Omth) call bounce_fast(v, eta, taub, bounceavg, timestep_transport, istate_dv) if (istate_dv == -1) then call error(fmt_dbg('VODE MXSTEP: mth=', dble(mth), ' ux=', ux, ' eta=', eta, ' taub=', taub)) else if (istate_dv /= 2) then call warning(fmt_dbg('dvode istate=', dble(istate_dv), ' at mth=', dble(mth), ' ux=', ux, ' eta=', eta)) else - if (abs(eta - etatp) < 1d-8*etatp) then + if (abs(eta - etatp) < 1.0e-8_dp*etatp) then call trace(fmt_dbg('near etatp: mth=', dble(mth), ' ux=', ux, ' eta=', eta, ' taub=', taub)) end if end if - Hmn2 = (bounceavg(3)**2 + bounceavg(4)**2)*(mi*(ux*vth)**2/2d0)**2 + Hmn2 = (bounceavg(3)**2 + bounceavg(4)**2)*(mi*(ux*vth)**2/2.0_dp)**2 attenuation_factor = nonlinear_attenuation(ux, eta, bounceavg, Omth, & dOmthdv, dOmthdeta, Hmn2) @@ -126,9 +127,9 @@ subroutine compute_transport_integral(vmin, vmax, vsteps, D, T) ux = ux + du end do - Dp = pi*vth**3/(16d0*R0*iota*(qi*B0/(mi*c))**2) - dsdreff = 2d0/a*sqrt(s) ! TODO: Use exact value instead of this approximation - D = dsdreff**(-2)*D/Dp + D_plateau = pi*vth**3/(16.0_dp*R0*iota*(qi*B0/(mi*c))**2) + dsdreff = 2.0_dp/a*sqrt(s) ! TODO: Use exact value instead of this approximation + D = dsdreff**(-2)*D/D_plateau call debug("compute_transport_integral complete") @@ -141,21 +142,21 @@ subroutine timestep_transport(v, eta, neq, t, y, ydot) ! More integrands may be added starting from y(3) ! - real(8), intent(in) :: v, eta + real(dp), intent(in) :: v, eta integer, intent(in) :: neq - real(8), intent(in) :: t - real(8), intent(in) :: y(neq) - real(8), intent(out) :: ydot(neq) + real(dp), intent(in) :: t + real(dp), intent(in) :: y(neq) + real(dp), intent(out) :: ydot(neq) ! BEGIN TODO: remove all of this after refactoring and re-use routine in orbit ! for y(1:3) - real(8) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3), Om_tB_v - real(8) :: t0 - complex(8) :: epsn, Hn ! relative amplitude of perturbation field epsn=Bn/B0 + real(dp) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3), Om_tB_v + real(dp) :: t0 + complex(dp) :: epsn, Hn ! relative amplitude of perturbation field epsn=Bn/B0 ! and Hamiltonian Hn = (H - H0)_n x(1) = s - x(2) = 0d0 + x(2) = 0.0_dp x(3) = y(1) call do_magfie(x, bmod, sqrtg, hder, hcovar, hctrvr, hcurl) call poloidal_velocity(v, eta, bmod, hctrvr(3), hder(3), y(2), ydot) @@ -170,20 +171,20 @@ subroutine timestep_transport(v, eta, neq, t, y, ydot) if (eta > etatp) then !t0 = 0.25*2*pi/Omth ! Different starting position in orbit - t0 = 0d0 - Hn = (2d0 - eta*bmod)*epsn*exp(imun*(q*mph*(y(1)) - mth*(t - t0)*Omth)) + t0 = 0.0_dp + Hn = (2.0_dp - eta*bmod)*epsn*exp(imun*(q*mph*(y(1)) - mth*(t - t0)*Omth)) else - Hn = (2d0 - eta*bmod)*epsn*exp(imun*(q*mph*(y(1)) - (mth + q*mph)*t*Omth)) + Hn = (2.0_dp - eta*bmod)*epsn*exp(imun*(q*mph*(y(1)) - (mth + q*mph)*t*Omth)) end if ydot(3) = real(Hn) ydot(4) = aimag(Hn) ! evaluate orbit averages for nonlinear attenuation if (nonlin) then - ydot(5) = 1d0/bmod + ydot(5) = 1.0_dp/bmod ydot(6) = bmod else - ydot(5:6) = 0d0 + ydot(5:6) = 0.0_dp end if end subroutine timestep_transport diff --git a/src/util.f90 b/src/util.f90 index a4fa0c4f..8403c0ce 100644 --- a/src/util.f90 +++ b/src/util.f90 @@ -3,22 +3,24 @@ ! Christopher Albert, 2015 module util + use iso_fortran_env, only: dp => real64 + implicit none save - complex(8), parameter :: imun = (0d0, 1d0) + complex(dp), parameter :: imun = (0.0_dp, 1.0_dp) - real(8), parameter, public :: pi = 4*atan(1d0) + real(dp), parameter, public :: pi = 4*atan(1.0_dp) - real(8), parameter, public :: & - qe = 4.803204d-10, & ! elementary charge - !me = 9.109382d-28, & ! electron mass, - mu = 1.660538d-24, & ! 1u - c = 2.997925d+10, & ! speed of light - !kb = 1.381649d-16, & ! Boltzmann constant - eV = 1.602176d-12 ! 1 electron volt + real(dp), parameter, public :: & + qe = 4.803204e-10_dp, & ! elementary charge + !me = 9.109382e-28_dp, & ! electron mass, + mu = 1.660538e-24_dp, & ! 1u + c = 2.997925e+10_dp, & ! speed of light + !kb = 1.381649e-16_dp, & ! Boltzmann constant + eV = 1.602176e-12_dp ! 1 electron volt - real(8), public :: qi = 1d0*qe, mi = 2.014d0*mu + real(dp), public :: qi = 1.0_dp*qe, mi = 2.014_dp*mu ! Flux-surface dependent (set from plasma data at each s) !$omp threadprivate (qi, mi) @@ -27,15 +29,15 @@ module util subroutine disp(str, val) character(*) :: str - real(8) :: val + real(dp) :: val write (*, '("'//str//'" ES16.9E2)') val end subroutine disp function linspace(a, b, cnt) integer :: cnt - real(8) :: linspace(cnt) + real(dp) :: linspace(cnt) integer :: i - real(8) :: a, b, delta + real(dp) :: a, b, delta delta = (b - a)/(cnt - 1) linspace = a + delta*(/(i, i=0, cnt - 1)/) @@ -69,8 +71,8 @@ end function newunit subroutine readdata(filename, ncol, data) character(*), intent(in) :: filename integer, intent(in) :: ncol - real(8), allocatable, intent(out) :: data(:, :) - real(8) :: buf(ncol) + real(dp), allocatable, intent(out) :: data(:, :) + real(dp) :: buf(ncol) integer :: lun, nrow, io, k open (unit=newunit(lun), file=filename, status="old") @@ -103,12 +105,12 @@ subroutine clearfile(filename) subroutine cached_spline(x, x_prev, spl_coeff, y) use spline, only: spline_val_0 - real(8), intent(in) :: x, x_prev - real(8), intent(in) :: spl_coeff(:, :, :) - real(8), intent(inout) :: y(:, :) + real(dp), intent(in) :: x, x_prev + real(dp), intent(in) :: spl_coeff(:, :, :) + real(dp), intent(inout) :: y(:, :) - real(8), parameter :: tol = 1.0d-12 - real(8) :: spl_val(3) + real(dp), parameter :: tol = 1.0e-12_dp + real(dp) :: spl_val(3) integer :: j if (abs(x - x_prev) < tol) return diff --git a/test/test_bounce.f90 b/test/test_bounce.f90 index 51796138..d69dd97e 100644 --- a/test/test_bounce.f90 +++ b/test/test_bounce.f90 @@ -1,10 +1,11 @@ program test_bounce_program + use iso_fortran_env, only: dp => real64 use util, only: qe, mu implicit none - real(8), parameter :: tol = 1e-13 - real(8) :: v, eta + real(dp), parameter :: tol = 1.0e-13_dp + real(dp) :: v, eta call setup call test_bounce @@ -30,7 +31,7 @@ subroutine setup sign_vpar = 1 v = vth - eta = 0.5d0*(etamin + etamax) + eta = 0.5_dp*(etamin + etamax) end subroutine setup @@ -39,14 +40,14 @@ subroutine setup_control mph, mth, magdrift, nopassing, pertfile, & nonlin, bfac, efac, inp_swi use neort_orbit, only: noshear - real(8) :: qs, ms - - s = 0.153d0 - M_t = 0.1d0 - qs = 1.0d0 - ms = 2.014d0 - vth = 37280978.0d0 - epsmn = 1e-3 + real(dp) :: qs, ms + + s = 0.153_dp + M_t = 0.1_dp + qs = 1.0_dp + ms = 2.014_dp + vth = 37280978.0_dp + epsmn = 1.0e-3_dp m0 = 0 mph = 18 mth = -1 @@ -55,8 +56,8 @@ subroutine setup_control noshear = .true. pertfile = .false. nonlin = .false. - bfac = 1.0d0 - efac = 1.0d0 + bfac = 1.0_dp + efac = 1.0_dp inp_swi = 8 M_t = M_t*efac/bfac @@ -67,12 +68,12 @@ end subroutine setup_control subroutine test_bounce use neort_orbit, only: bounce, nvar - real(8) :: taub, bounceavg(nvar), bounceavg_tmp(nvar) + real(dp) :: taub, bounceavg(nvar), bounceavg_tmp(nvar) - bounceavg = 0d0 + bounceavg = 0.0_dp call bounce(v, eta, taub, bounceavg_tmp) - bounceavg = 0d0 + bounceavg = 0.0_dp call bounce(v, eta, taub, bounceavg) if (maxval(abs(bounceavg - bounceavg_tmp)) > tol) then @@ -86,18 +87,18 @@ end subroutine test_bounce subroutine test_bounce_fast use neort_orbit, only: bounce, bounce_fast, nvar, timestep - real(8) :: taub, bounceavg(nvar), bounceavg_tmp(nvar) + real(dp) :: taub, bounceavg(nvar), bounceavg_tmp(nvar) - bounceavg = 0d0 + bounceavg = 0.0_dp - bounceavg = 0d0 + bounceavg = 0.0_dp call bounce(v, eta, taub, bounceavg_tmp, taub) - bounceavg = 0d0 + bounceavg = 0.0_dp call bounce_fast(v, eta, taub, bounceavg, timestep) - if (maxval(abs((bounceavg - bounceavg_tmp)/(bounceavg+1d-6))) & - > 1d-3) then + if (maxval(abs((bounceavg - bounceavg_tmp)/(bounceavg+1.0e-6_dp))) & + > 1.0e-3_dp) then print *, 'test_bounce_fast failed' error stop end if diff --git a/test/test_frequencies.f90 b/test/test_frequencies.f90 index 611828e8..f727575b 100644 --- a/test/test_frequencies.f90 +++ b/test/test_frequencies.f90 @@ -1,4 +1,5 @@ program test_frequencies + use iso_fortran_env, only: dp => real64 use util use do_magfie_mod, only: s, R0, a, inp_swi, do_magfie_init, do_magfie use neort, only: init @@ -9,14 +10,14 @@ program test_frequencies implicit none integer, parameter :: neta = 1000 ! Steps in eta - real(8), parameter :: M_t = 8.2512439552d-02 ! Unscaled Mach number + real(dp), parameter :: M_t = 8.2512439552e-02_dp ! Unscaled Mach number integer :: i, fid - real(8) :: Omth, dOmthdv, dOmthdeta, Omph, dOmphdv, dOmphdeta - real(8), parameter :: scalfac_energy = 1.0d0 - real(8), parameter :: scalfac_efield = 1.0d0 + real(dp) :: Omth, dOmthdv, dOmthdeta, Omph, dOmphdv, dOmphdeta + real(dp), parameter :: scalfac_energy = 1.0_dp + real(dp), parameter :: scalfac_efield = 1.0_dp - real(8) :: v, eta + real(dp) :: v, eta inp_swi = 9 ! ASDEX Upgrade format s = 0.17 ! Normalized toroidal flux @@ -27,7 +28,7 @@ program test_frequencies mth = 1 ! Canonical poloidal harmonic mph = 2 ! Canonical toroidal harmonic - vth = 4.5972926215d7 ! Unscaled thermal velocity + vth = 4.5972926215e7_dp ! Unscaled thermal velocity Om_tE = vth*M_t/(R0*scalfac_efield) ! ExB precession frequency vth = vth/sqrt(scalfac_energy) @@ -37,7 +38,7 @@ program test_frequencies open(newunit=fid, file='canonical_freqs_vs_eta_t.dat') write(fid,*) '%# eta [1/G], omega_b [rad/s], Omega_tor [rad/s] trapped' - sign_vpar = 1d0 + sign_vpar = 1.0_dp do i = 1,neta eta = etatp + i*(etadt-etatp)/neta call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) @@ -48,7 +49,7 @@ program test_frequencies open(newunit=fid, file='canonical_freqs_vs_eta_pco.dat') write(fid,*) '%# eta [1/G], omega_b [rad/s], Omega_tor [rad/s] co-passing' - sign_vpar = 1d0 + sign_vpar = 1.0_dp do i = 0,neta-1 eta = i*etatp/neta call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) @@ -59,7 +60,7 @@ program test_frequencies open(newunit=fid, file='canonical_freqs_vs_eta_pct.dat') write(fid,*) '%# eta [1/G], omega_b [rad/s], Omega_tor [rad/s] co-passing' - sign_vpar = -1d0 + sign_vpar = -1.0_dp do i = 0,neta-1 eta = i*etatp/neta call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) diff --git a/test/test_misc.f90 b/test/test_misc.f90 index 441c0cb2..cde8482b 100644 --- a/test/test_misc.f90 +++ b/test/test_misc.f90 @@ -1,4 +1,5 @@ program test_misc + use iso_fortran_env, only: dp => real64 use util, only: readdata, disp, c, qi, mi, pi use driftorbit, only: etamin, etamax, etatp, etadt, epsst_spl, epst_spl, epssp_spl, & epsp_spl, sign_vpar, mth, nlev, nopassing, epsp, epst, dVds @@ -24,28 +25,28 @@ program test_misc subroutine test_torfreq integer, parameter :: n = 100 integer :: k - real(8) :: Omph, dOmphdv, dOmphdeta, dOmphds - real(8) :: Omth, dOmthdv, dOmthdeta, dOmthds - real(8) :: OmtB, dOmtBdv, dOmtBdeta - real(8) :: aa, b - real(8) :: taub, bounceavg(nvar) - real(8) :: v, eta + real(dp) :: Omph, dOmphdv, dOmphdeta, dOmphds + real(dp) :: Omth, dOmthdv, dOmthdeta, dOmthds + real(dp) :: OmtB, dOmtBdv, dOmtBdeta + real(dp) :: aa, b + real(dp) :: taub, bounceavg(nvar) + real(dp) :: v, eta v = vth etamin = etatp - etamax = etatp + (etadt - etatp)*(1d0 - epsst_spl) + etamax = etatp + (etadt - etatp)*(1.0_dp - epsst_spl) b = log(epst_spl) - aa = 1d0/(n - 1d0)*(log(etamax/etamin - 1d0) - b) + aa = 1.0_dp/(n - 1.0_dp)*(log(etamax/etamin - 1.0_dp) - b) eta = etamax sign_vpar = 1 call disp("test_torfreq: vth = ", vth) call disp("test_torfreq: v/vth = ", v/vth) - call disp("test_torfreq: mph = ", 1d0*mph) - call disp("test_torfreq: mth = ", 1d0*mth) + call disp("test_torfreq: mph = ", 1.0_dp*mph) + call disp("test_torfreq: mth = ", 1.0_dp*mth) call disp("test_torfreq: Om_tE = ", Om_tE) call disp("test_torfreq: Om_tB_ref = ", c*mi*vth**2/(2*qi*psi_pr)) call bounce(v, eta, taub, bounceavg) @@ -54,7 +55,7 @@ subroutine test_torfreq call disp("test_torfreq: etamax = ", etamax) call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) - call disp("test_torfreq: Om_th_approx = ", v/(q*R0*sqrt(2d0/eps))) + call disp("test_torfreq: Om_th_approx = ", v/(q*R0*sqrt(2.0_dp/eps))) call disp("test_torfreq: Om_th_deeptrap = ", Omth) open (unit=9, file=trim(adjustl(runname))//"_torfreq.out", recl=1024) @@ -74,11 +75,11 @@ subroutine test_torfreq "14:dOmphdeta "// & "15:dOmphds " do k = 0, n - 1 - eta = etamin*(1d0 + exp(aa*k + b)) + eta = etamin*(1.0_dp + exp(aa*k + b)) call Om_ph(v, eta, Omph, dOmphdv, dOmphdeta) call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) call Om_tB(v, eta, OmtB, dOmtBdv, dOmtBdeta) - call d_Om_ds(v, eta, 2d0*pi/Omth, dOmthds, dOmphds) + call d_Om_ds(v, eta, 2.0_dp*pi/Omth, dOmthds, dOmphds) write (9, *) eta, etatp, etadt, & Om_tE, OmtB, dOmtbdv, dOmtbdeta, & Omth, dOmthdv, dOmthdeta, dOmthds, & @@ -90,15 +91,15 @@ end subroutine test_torfreq subroutine test_resline integer, parameter :: n = 500 integer :: k - real(8) :: vmin, vmax, v - real(8) :: etarest(2), etaresp(2) - real(8) :: roots(nlev, 3) + real(dp) :: vmin, vmax, v + real(dp) :: etarest(2), etaresp(2) + real(dp) :: roots(nlev, 3) integer :: nroots, kr print *, "test_resline" - vmin = 1d-6*vth - vmax = 10d0*vth + vmin = 1.0e-6_dp*vth + vmax = 10.0_dp*vth etaresp = etatp etarest = etatp @@ -108,21 +109,21 @@ subroutine test_resline open (unit=10, file=trim(adjustl(runname))//"_resline_pct.out", recl=1024) open (unit=11, file=trim(adjustl(runname))//"_resline_t.out", recl=1024) do k = 0, n - 1 - v = vmin + k/(n - 1d0)*(vmax - vmin) + v = vmin + k/(n - 1.0_dp)*(vmax - vmin) if (.not. nopassing) then ! resonance (passing) sign_vpar = 1 call driftorbit_coarse(v, etatp*epsp, etatp*(1 - epsp), roots, nroots) do kr = 1, nroots - etaresp = driftorbit_root(v, 1d-8*abs(Om_tE), roots(kr, 1), roots(kr, 2)) - write (9, *) v/vth, kr, etaresp(1), 0d0, etatp + etaresp = driftorbit_root(v, 1.0e-8_dp*abs(Om_tE), roots(kr, 1), roots(kr, 2)) + write (9, *) v/vth, kr, etaresp(1), 0.0_dp, etatp end do sign_vpar = -1 call driftorbit_coarse(v, etatp*epsp, etatp*(1 - epsp), roots, nroots) do kr = 1, nroots - etaresp = driftorbit_root(v, 1d-8*abs(Om_tE), roots(kr, 1), roots(kr, 2)) - write (10, *) v/vth, kr, etaresp(1), 0d0, etatp + etaresp = driftorbit_root(v, 1.0e-8_dp*abs(Om_tE), roots(kr, 1), roots(kr, 2)) + write (10, *) v/vth, kr, etaresp(1), 0.0_dp, etatp end do end if @@ -130,7 +131,7 @@ subroutine test_resline sign_vpar = 1 call driftorbit_coarse(v, etatp*(1 + epst), etadt*(1 - epst), roots, nroots) do kr = 1, nroots - etarest = driftorbit_root(v, 1d-8*abs(Om_tE), roots(kr, 1), roots(kr, 2)) + etarest = driftorbit_root(v, 1.0e-8_dp*abs(Om_tE), roots(kr, 1), roots(kr, 2)) write (11, *) v/vth, kr, etarest(1), etatp, etadt end do end do @@ -141,7 +142,7 @@ end subroutine test_resline subroutine test_profile integer :: k - real(8), allocatable :: data(:, :) + real(dp), allocatable :: data(:, :) print *, "test_profile" @@ -153,15 +154,15 @@ subroutine test_profile s = data(k, 1) call output_flux_surface_data(9) end do - s = 0.96 + s = 0.96_dp call output_flux_surface_data(9) - s = 0.97 + s = 0.97_dp call output_flux_surface_data(9) - s = 0.98 + s = 0.98_dp call output_flux_surface_data(9) - s = 0.99 + s = 0.99_dp call output_flux_surface_data(9) - s = 1.0 + s = 1.0_dp call output_flux_surface_data(9) close (unit=9) @@ -171,10 +172,10 @@ end subroutine test_profile subroutine output_flux_surface_data(unit) integer, intent(in) :: unit - real(8) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) + real(dp) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) x(1) = s - x(2) = 0d0 - x(3) = 0d0 + x(2) = 0.0_dp + x(3) = 0.0_dp call read_and_init_plasma_input("plasma.in", s) call init_flux_surface_average(s) call do_magfie(x, bmod, sqrtg, hder, hcovar, hctrvr, hcurl) diff --git a/test/test_neort_lib.f90 b/test/test_neort_lib.f90 index 34773b62..5a90eec3 100644 --- a/test/test_neort_lib.f90 +++ b/test/test_neort_lib.f90 @@ -20,7 +20,7 @@ program test_neort_lib subroutine test_single_s() type(transport_data_t) :: result - real(dp), parameter :: TEST_S = 0.5d0 + real(dp), parameter :: TEST_S = 0.5_dp print *, "Test: single s computation..." @@ -28,7 +28,7 @@ subroutine test_single_s() call neort_prepare_splines("plasma.in", "profile.in") call neort_compute_at_s(TEST_S, result) - if (abs(result%torque%s - TEST_S) > 1d-12) then + if (abs(result%torque%s - TEST_S) > 1.0e-12_dp) then print *, "ERROR: s value mismatch, expected", TEST_S, "got", result%torque%s error stop "test_single_s failed" end if @@ -37,7 +37,7 @@ subroutine test_single_s() error stop "test_single_s: expected torque computation" end if - if (result%torque%dVds <= 0d0) then + if (result%torque%dVds <= 0.0_dp) then error stop "test_single_s: dVds should be positive" end if @@ -49,7 +49,7 @@ subroutine test_parallel_loop() integer, parameter :: N_THREADS = 4 type(transport_data_t) :: results(N_S) real(dp) :: s_array(N_S) - real(dp), parameter :: TOL = 1d-10 + real(dp), parameter :: TOL = 1.0e-10_dp integer :: i logical :: all_positive @@ -58,7 +58,7 @@ subroutine test_parallel_loop() call omp_set_num_threads(N_THREADS) do i = 1, N_S - s_array(i) = 0.2d0 + 0.1d0 * real(i, dp) + s_array(i) = 0.2_dp + 0.1_dp * real(i, dp) end do call neort_init("driftorbit", "in_file") @@ -76,7 +76,7 @@ subroutine test_parallel_loop() print *, "ERROR: s mismatch at index", i error stop "test_parallel_loop: s value mismatch" end if - if (results(i)%torque%dVds <= 0d0) then + if (results(i)%torque%dVds <= 0.0_dp) then all_positive = .false. end if end do @@ -93,7 +93,7 @@ subroutine test_spline_update() use util, only: readdata type(transport_data_t) :: result1, result2 - real(dp), parameter :: TEST_S = 0.5d0 + real(dp), parameter :: TEST_S = 0.5_dp integer :: nplasma real(dp) :: am1, am2, Z1, Z2 real(dp), allocatable :: plasma_data(:, :), profile_data(:, :) @@ -111,11 +111,11 @@ subroutine test_spline_update() call neort_prepare_splines(nplasma, am1, am2, Z1, Z2, plasma_data, profile_data) call neort_compute_at_s(TEST_S, result2) - if (abs(result1%torque%s - result2%torque%s) > 1d-12) then + if (abs(result1%torque%s - result2%torque%s) > 1.0e-12_dp) then error stop "test_spline_update: results should match" end if - if (abs(result1%torque%Tt - result2%torque%Tt) > 1d-10) then + if (abs(result1%torque%Tt - result2%torque%Tt) > 1.0e-10_dp) then error stop "test_spline_update: torque should match" end if diff --git a/test/test_omega_prime.f90 b/test/test_omega_prime.f90 index d4610de8..55fd29bc 100644 --- a/test/test_omega_prime.f90 +++ b/test/test_omega_prime.f90 @@ -1,4 +1,5 @@ program test_omage_prime_prog + use iso_fortran_env, only: dp => real64 use do_magfie_mod, only: do_magfie_init, params0, sign_theta, s use do_magfie_pert_mod, only: do_magfie_pert_init use neort, only: check_magfie, write_magfie_data_to_files, init_profiles, & @@ -16,27 +17,27 @@ program test_omage_prime_prog use driftorbit implicit none - real(8) :: DELTA = 1d-8 + real(dp) :: DELTA = 1.0e-8_dp type :: freq_data_t - real(8) :: s, v, eta - real(8) :: J(3), Jbar(3) - real(8) :: Om, Ompr_old, Ompr_new - real(8) :: dOmds, dOmdv, dOmdeta, dOmdpph - real(8) :: Omth, dOmthds, dOmthdv, dOmthdeta - real(8) :: Omph, dOmphds, dOmphdv, dOmphdeta - real(8) :: Om_tE, dOm_tEds + real(dp) :: s, v, eta + real(dp) :: J(3), Jbar(3) + real(dp) :: Om, Ompr_old, Ompr_new + real(dp) :: dOmds, dOmdv, dOmdeta, dOmdpph + real(dp) :: Omth, dOmthds, dOmthdv, dOmthdeta + real(dp) :: Omph, dOmphds, dOmphdv, dOmphdeta + real(dp) :: Om_tE, dOm_tEds end type - real(8), allocatable :: spl_psi_pol(:, :) + real(dp), allocatable :: spl_psi_pol(:, :) call main contains subroutine main integer, parameter :: NUM_SAMPLES = 32 - real(8) :: s0, v0, eta0 - real(8) :: eta_res(2) + real(dp) :: s0, v0, eta0 + real(dp) :: eta_res(2) type(freq_data_t) :: freq_data(NUM_SAMPLES) @@ -44,10 +45,10 @@ subroutine main call read_and_set_config(runname) call setup s0 = s - v0 = 0.9d0*vth + v0 = 0.9_dp*vth eta_res = first_resonance(v0) eta0 = eta_res(1) - !eta0 = 0.5d0*(etamin + etamax) + !eta0 = 0.5_dp*(etamin + etamax) call spline_psi_pol call sample_values(s0, v0, eta0, freq_data) @@ -55,11 +56,11 @@ subroutine main end subroutine main subroutine sample_values(s0, v0, eta0, freq_data) - real(8), intent(in) :: s0, v0, eta0 + real(dp), intent(in) :: s0, v0, eta0 type(freq_data_t), intent(out) :: freq_data(:) integer :: i - real(8) :: xi + real(dp) :: xi freq_data(1) % s = s0 s = freq_data(1) % s @@ -68,12 +69,12 @@ subroutine sample_values(s0, v0, eta0, freq_data) call test_omega_prime(freq_data(1)) do i = 2, size(freq_data) call random_number(xi) - freq_data(i) % s = s0*(1d0 + DELTA*(xi - 0.5d0)) + freq_data(i) % s = s0*(1.0_dp + DELTA*(xi - 0.5_dp)) s = freq_data(i) % s call random_number(xi) - freq_data(i) % v = v0*(1d0 + DELTA*(xi - 0.5d0)) + freq_data(i) % v = v0*(1.0_dp + DELTA*(xi - 0.5_dp)) call random_number(xi) - freq_data(i) % eta = eta0 + (etamax - etamin)*DELTA*(xi - 0.5d0) + freq_data(i) % eta = eta0 + (etamax - etamin)*DELTA*(xi - 0.5_dp) call test_omega_prime(freq_data(i)) end do end subroutine sample_values @@ -118,17 +119,17 @@ subroutine setup end subroutine setup function first_resonance(v) result(eta_res) - real(8), intent(in) :: v - real(8) :: eta_res(2) - real(8) :: roots(nlev, 3) + real(dp), intent(in) :: v + real(dp) :: eta_res(2) + real(dp) :: roots(nlev, 3) integer :: nroots call driftorbit_coarse(v, etamin, etamax, roots, nroots) - eta_res = driftorbit_root(v, 1d-8*abs(Om_tE), roots(1, 1), roots(1, 2)) + eta_res = driftorbit_root(v, 1.0e-8_dp*abs(Om_tE), roots(1, 1), roots(1, 2)) end function first_resonance subroutine test_omega_prime(f) type(freq_data_t), intent(inout) :: f - real(8) :: bounceavg(nvar) + real(dp) :: bounceavg(nvar) call setup @@ -136,7 +137,7 @@ subroutine test_omega_prime(f) print *, "dOmthds_test: ", dOmthds_test(f % v, f % eta) - call bounce_fast(f % v, f % eta, 2d0*pi/abs(f % Omth), bounceavg, timestep_transport) + call bounce_fast(f % v, f % eta, 2.0_dp*pi/abs(f % Omth), bounceavg, timestep_transport) f % Ompr_old = omega_prime_old(f % v/vth, f % eta, f % Omth, f % dOmdv, f % dOmdeta, f % dOmdpph) f % Ompr_new = omega_prime(f % v/vth, f % eta, bounceavg, f % Omth, f % dOmdv, f % dOmdeta, f % dOmdpph) @@ -152,10 +153,10 @@ end subroutine test_omega_prime subroutine compute_frequencies(f) type(freq_data_t), intent(inout) :: f - real(8) :: taub, bounceavg(nvar) + real(dp) :: taub, bounceavg(nvar) call Om_th(f % v, f % eta, f % Omth, f % dOmthdv, f % dOmthdeta) - taub = 2d0*pi/abs(f % Omth) + taub = 2.0_dp*pi/abs(f % Omth) call bounce_fast(f % v, f % eta, taub, bounceavg, timestep_transport) call Om_ph(f % v, f % eta, f % Omph, f % dOmphdv, f % dOmphdeta) call d_Om_ds(f % v, f % eta, taub, f % dOmthds, f % dOmphds) @@ -169,24 +170,24 @@ subroutine compute_frequencies(f) end subroutine compute_frequencies subroutine compute_invariants(v, eta, J) - real(8), intent(in) :: v, eta - real(8), intent(out) :: J(3) + real(dp), intent(in) :: v, eta + real(dp), intent(out) :: J(3) integer, parameter :: neq = 3 - real(8) :: taub, bounceavg(neq) - real(8) :: bounceint(neq + 1) - real(8) :: bmod, htheta - real(8) :: y0(neq) + real(dp) :: taub, bounceavg(neq) + real(dp) :: bounceint(neq + 1) + real(dp) :: bmod, htheta + real(dp) :: y0(neq) ! Initialize bounce-averated quantities y0. Their meaning ! is defined inside subroutine timestep (thin orbit integration) call evaluate_bfield_local(bmod, htheta) y0(1) = th0 - y0(2) = sign(1d0, htheta)*sign_vpar*vpar(v, eta, bmod) - y0(3) = 1d-15 + y0(2) = sign(1.0_dp, htheta)*sign_vpar*vpar(v, eta, bmod) + y0(3) = 1.0e-15_dp - taub = 2.0*pi/(vperp(v, eta, bmod)*iota/R0*sqrt(eps/2d0)) - bounceint = bounce_integral(v, eta, neq, y0, taub/5d0, timestep_invariants) + taub = 2.0_dp*pi/(vperp(v, eta, bmod)*iota/R0*sqrt(eps/2.0_dp)) + bounceint = bounce_integral(v, eta, neq, y0, taub/5.0_dp, timestep_invariants) taub = bounceint(1) bounceavg = bounceint(2:)/taub @@ -204,15 +205,15 @@ subroutine compute_invariants(v, eta, J) end subroutine compute_invariants subroutine timestep_invariants(v, eta, neq, t, y, ydot) - real(8), intent(in) :: v, eta, t + real(dp), intent(in) :: v, eta, t integer, intent(in) :: neq - real(8), intent(in) :: y(neq) - real(8), intent(out) :: ydot(neq) + real(dp), intent(in) :: y(neq) + real(dp), intent(out) :: ydot(neq) - real(8) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) + real(dp) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) x(1) = s - x(2) = 0d0 + x(2) = 0.0_dp x(3) = y(1) call do_magfie(x, bmod, sqrtg, hder, hcovar, hctrvr, hcurl) @@ -222,30 +223,30 @@ subroutine timestep_invariants(v, eta, neq, t, y, ydot) end subroutine timestep_invariants pure function Jperp(v, eta) - real(8) :: Jperp - real(8), intent(in) :: v, eta - Jperp = 0.5d0*mi*v**2*mi*c/qi*eta + real(dp) :: Jperp + real(dp), intent(in) :: v, eta + Jperp = 0.5_dp*mi*v**2*mi*c/qi*eta end function Jperp function pphi() - real(8) :: pphi + real(dp) :: pphi pphi = -qi/c*psi_pol() end function pphi function psi_pol() - real(8) :: psi_pol, psi_pol_result(3) + real(dp) :: psi_pol, psi_pol_result(3) psi_pol_result = spline_val_0(spl_psi_pol, s) psi_pol = psi_pol_result(1) end function psi_pol subroutine psi_pol_grid(psi) - real(8), allocatable, intent(out) :: psi(:) + real(dp), allocatable, intent(out) :: psi(:) integer :: i associate (sx => params0(:, 1), iota => params0(:, 2)) allocate (psi(size(sx))) - psi(1) = 0d0 + psi(1) = 0.0_dp do i = 2, size(sx) psi(i) = psi(i - 1) & + psi_pr*trapz_step(sx(i - 1), sx(i), iota(i - 1), iota(i)) @@ -254,7 +255,7 @@ subroutine psi_pol_grid(psi) end subroutine psi_pol_grid subroutine spline_psi_pol - real(8), allocatable :: psi(:) + real(dp), allocatable :: psi(:) associate (sx => params0(:, 1)) call psi_pol_grid(psi) @@ -264,48 +265,48 @@ subroutine spline_psi_pol end subroutine spline_psi_pol function trapz(x, y) - real(8), intent(in) :: x(:), y(:) - real(8) :: trapz + real(dp), intent(in) :: x(:), y(:) + real(dp) :: trapz integer :: i - trapz = 0d0 + trapz = 0.0_dp do i = 2, size(x) trapz = trapz + trapz_step(x(i - 1), x(i), y(i - 1), y(i)) end do end function trapz function trapz_step(x1, x2, y1, y2) - real(8), intent(in) :: x1, x2, y1, y2 - real(8) :: trapz_step + real(dp), intent(in) :: x1, x2, y1, y2 + real(dp) :: trapz_step - trapz_step = 0.5d0*(x2 - x1)*(y1 + y2) + trapz_step = 0.5_dp*(x2 - x1)*(y1 + y2) end function trapz_step function dOmthds_test(v, eta) - real(8) :: dOmthds_test - real(8), intent(in) :: v, eta + real(dp) :: dOmthds_test + real(dp), intent(in) :: v, eta - real(8) :: s0, ds - real(8) :: Omth1, Omth2 + real(dp) :: s0, ds + real(dp) :: Omth1, Omth2 ds = DELTA s0 = s s = s0 + ds - Omth1 = 2d0*pi/bounce_time(v, eta) + Omth1 = 2.0_dp*pi/bounce_time(v, eta) s = s0 - ds - Omth2 = 2d0*pi/bounce_time(v, eta) + Omth2 = 2.0_dp*pi/bounce_time(v, eta) - dOmthds_test = (Omth1 - Omth2)/(2d0*ds) + dOmthds_test = (Omth1 - Omth2)/(2.0_dp*ds) s = s0 end function dOmthds_test function omega_prime_old(ux, eta, Omth, dOmdv, dOmdeta, dOmdpph) - real(8), intent(in) :: ux, eta, Omth, dOmdv, dOmdeta, dOmdpph - real(8) :: omega_prime_old - omega_prime_old = mth*(eta*dOmdeta - ux*vth/2*dOmdv)/(mi*(ux*vth)**2/(2d0*Omth)) + dOmdpph + real(dp), intent(in) :: ux, eta, Omth, dOmdv, dOmdeta, dOmdpph + real(dp) :: omega_prime_old + omega_prime_old = mth*(eta*dOmdeta - ux*vth/2*dOmdv)/(mi*(ux*vth)**2/(2.0_dp*Omth)) + dOmdpph end function omega_prime_old end program test_omage_prime_prog diff --git a/test/test_parallel.f90 b/test/test_parallel.f90 index cfcb56cc..e025efa8 100644 --- a/test/test_parallel.f90 +++ b/test/test_parallel.f90 @@ -18,8 +18,8 @@ program test_parallel implicit none integer, parameter :: N_THREADS = 8 - real(dp), parameter :: TOL_INIT = 1.0d-12 - real(dp), parameter :: TOL_TRANSPORT = 1.0d-10 + real(dp), parameter :: TOL_INIT = 1.0e-12_dp + real(dp), parameter :: TOL_TRANSPORT = 1.0e-10_dp type :: init_results_t real(dp) :: s, psi_pr, Bthcov, Bphcov, q, iota, eps, B0h @@ -63,7 +63,7 @@ subroutine run_test() print *, "Input files read." print * - test_s = 0.5d0 + test_s = 0.5_dp !$omp parallel private(thread_id, x, bmod, sqrtg, bder, hcovar, hctrvr, hcurl) thread_id = omp_get_thread_num() + 1 @@ -208,7 +208,7 @@ subroutine check_value(name, expected, actual, thread_num, tol, all_match) real(dp) :: diff, rel_diff diff = abs(expected - actual) - if (abs(expected) > 1d-30) then + if (abs(expected) > 1.0e-30_dp) then rel_diff = diff / abs(expected) else rel_diff = diff diff --git a/test/test_reslines.f90 b/test/test_reslines.f90 index 2bb67d68..18151cdc 100644 --- a/test/test_reslines.f90 +++ b/test/test_reslines.f90 @@ -1,4 +1,5 @@ program test_reslines + use iso_fortran_env, only: dp => real64 use util use do_magfie_mod, only: s, R0, a, inp_swi, do_magfie_init, do_magfie use do_magfie_pert_mod, only: mph @@ -9,13 +10,13 @@ program test_reslines use neort_freq, only: Om_th, Om_ph implicit none - real(8), parameter :: vnorm = sqrt(0.15) + real(dp), parameter :: vnorm = sqrt(0.15d0) integer :: fid integer, parameter :: ns = 250 integer :: ks - real(8) :: profile(ns, 3) + real(dp) :: profile(ns, 3) integer, parameter :: nm2 = 11, nm3 = 2 integer :: km2, km3 @@ -56,40 +57,40 @@ program test_reslines subroutine resline - real(8) :: v, eta_res(2) - real(8) :: roots(nlev, 3) + real(dp) :: v, eta_res(2) + real(dp) :: roots(nlev, 3) integer :: nroots, kr v = vnorm*vth ! Trapped - sign_vpar = 1d0 + sign_vpar = 1.0_dp etamin = (1+epst)*etatp etamax = (1-epst)*etadt call driftorbit_coarse(v, etamin, etamax, roots, nroots) if(nroots == 0) return do kr = 1,nroots - eta_res = driftorbit_root(v, 1d-8*abs(Om_tE), roots(kr,1), roots(kr,2)) + eta_res = driftorbit_root(v, 1.0e-8_dp*abs(Om_tE), roots(kr,1), roots(kr,2)) write(fid,*) s, v/vth, eta_res(1), mth, mph, kr, 0 end do ! Co-passing - sign_vpar = 1d0 + sign_vpar = 1.0_dp etamin = epsp*etatp etamax = (1-epsp)*etatp call driftorbit_coarse(v, etamin, etamax, roots, nroots) if(nroots == 0) return do kr = 1,nroots - eta_res = driftorbit_root(v, 1d-8*abs(Om_tE), roots(kr,1), roots(kr,2)) + eta_res = driftorbit_root(v, 1.0e-8_dp*abs(Om_tE), roots(kr,1), roots(kr,2)) write(fid,*) s, v/vth, eta_res(1), mth, mph, kr, 1 end do ! Ctr-passing - sign_vpar = -1d0 + sign_vpar = -1.0_dp etamin = epsp*etatp etamax = (1-epsp)*etatp call driftorbit_coarse(v, etamin, etamax, roots, nroots) if(nroots == 0) return do kr = 1,nroots - eta_res = driftorbit_root(v, 1d-8*abs(Om_tE), roots(kr,1), roots(kr,2)) + eta_res = driftorbit_root(v, 1.0e-8_dp*abs(Om_tE), roots(kr,1), roots(kr,2)) write(fid,*) s, v/vth, eta_res(1), mth, mph, kr, -1 end do diff --git a/test/test_torque.f90 b/test/test_torque.f90 index d852e08a..e28bd23e 100644 --- a/test/test_torque.f90 +++ b/test/test_torque.f90 @@ -1,4 +1,5 @@ program test_torque_prog + use iso_fortran_env, only: dp => real64 use neort, only: check_magfie, init_profiles, & init, compute_transport_harmonic, write_magfie_data_to_files use neort_config, only: read_and_set_config @@ -10,14 +11,14 @@ program test_torque_prog use util implicit none - real(8) :: Dco(2), Dctr(2), Dt(2) - real(8) :: Tco, Tctr, Tt + real(dp) :: Dco(2), Dctr(2), Dt(2) + real(dp) :: Tco, Tctr, Tt integer, parameter :: MTH = -3 - real(8), parameter :: AVG_NABLA_S = 0.02162413513580247d0 - real(8), parameter :: DV_OVER_DS = 12621249.179676009d0 - real(8), parameter :: FLUX_SURFACE_AREA = DV_OVER_DS*AVG_NABLA_S - real(8), parameter :: TOL = 1d-7 + real(dp), parameter :: AVG_NABLA_S = 0.02162413513580247_dp + real(dp), parameter :: DV_OVER_DS = 12621249.179676009_dp + real(dp), parameter :: FLUX_SURFACE_AREA = DV_OVER_DS*AVG_NABLA_S + real(dp), parameter :: TOL = 1.0e-7_dp character(1024) :: TEST_RUN = "driftorbit64.new" type(magfie_data_t) :: magfie_data type(transport_harmonic_t) :: harmonic_data @@ -38,12 +39,12 @@ subroutine setup call check_magfie(magfie_data) call write_magfie_data_to_files(magfie_data, runname) - Dco = 0d0 - Dctr = 0d0 - Dt = 0d0 - Tco = 0d0 - Tctr = 0d0 - Tt = 0d0 + Dco = 0.0_dp + Dctr = 0.0_dp + Dt = 0.0_dp + Tco = 0.0_dp + Tctr = 0.0_dp + Tt = 0.0_dp call compute_transport_harmonic(MTH, Dco, Dctr, Dt, Tco, Tctr, Tt, harmonic_data) call correct_transport_coeff(Dco) call correct_transport_coeff(Dctr) @@ -54,7 +55,7 @@ subroutine setup end subroutine setup subroutine test_torque - real(8) :: TTco, TTctr, TTt + real(dp) :: TTco, TTctr, TTt TTco = torque_from_transport(Dco) TTctr = torque_from_transport(Dctr) @@ -68,11 +69,11 @@ subroutine test_torque print *, "TTctr/Tctr = ", TTctr/Tctr print *, "TTt/Tt = ", TTt/Tt - if (abs(TTctr/Tctr - 1d0) > TOL) then + if (abs(TTctr/Tctr - 1.0_dp) > TOL) then error stop end if - if (abs(TTt/Tt - 1d0) > TOL) then + if (abs(TTt/Tt - 1.0_dp) > TOL) then error stop end if @@ -81,35 +82,35 @@ subroutine test_torque end subroutine test_torque function torque_from_transport(D) result(Tphi) - real(8), intent(in) :: D(2) - real(8) :: Tphi - real(8) :: Dp - real(8) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) - real(8) :: sqrtgBth + real(dp), intent(in) :: D(2) + real(dp) :: Tphi + real(dp) :: D_plateau + real(dp) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3) + real(dp) :: sqrtgBth x(1) = s - x(2) = 0d0 - x(3) = 0d0 + x(2) = 0.0_dp + x(3) = 0.0_dp call do_magfie(x, bmod, sqrtg, hder, hcovar, hctrvr, hcurl) sqrtgBth = sqrtg*hctrvr(3)*bmod*AVG_NABLA_S print *, "sqrtgBth = ", sqrtgBth - Dp = pi*vth**3/(16d0*R0*iota*(qi*B0/(mi*c))**2) - Tphi = (sqrtgBth/c)*qe*(-ni1*Dp*(D(1)*A1 + D(2)*A2))*AVG_NABLA_S + D_plateau = pi*vth**3/(16.0_dp*R0*iota*(qi*B0/(mi*c))**2) + Tphi = (sqrtgBth/c)*qe*(-ni1*D_plateau*(D(1)*A1 + D(2)*A2))*AVG_NABLA_S Tphi = Tphi end function torque_from_transport ! Correct ad-hoc approximation for AVG_NABLA_S ! This then finally cancels out in torque calculation subroutine correct_transport_coeff(D) - real(8), intent(inout) :: D(2) - D = D*(2d0/a*sqrt(s))**2/AVG_NABLA_S**2 + real(dp), intent(inout) :: D(2) + D = D*(2.0_dp/a*sqrt(s))**2/AVG_NABLA_S**2 end subroutine correct_transport_coeff ! Correct torque calculation subroutine correct_torque(T) - real(8), intent(inout) :: T + real(dp), intent(inout) :: T T = T/DV_OVER_DS end subroutine correct_torque end program test_torque_prog From 78440c768826df890862620c364bb97068f5cdf3 Mon Sep 17 00:00:00 2001 From: Andreas Zach Date: Sat, 10 Jan 2026 01:10:13 +0100 Subject: [PATCH 2/3] Reformat touched lines --- src/attenuation_factor.f90 | 2 +- src/collis_nbi.f90 | 80 ++++++++++++++-------------- src/do_magfie_standalone.f90 | 100 +++++++++++++++++------------------ src/freq.f90 | 81 ++++++++++++++-------------- src/magfie.f90 | 10 ++-- src/neort.f90 | 16 +++--- src/nonlin.f90 | 36 ++++++------- src/orbit.f90 | 59 ++++++++++----------- src/profiles.f90 | 19 +++---- src/resonance.f90 | 20 +++---- src/transport.f90 | 59 +++++++++++---------- src/util.f90 | 12 ++--- test/test_bounce.f90 | 10 ++-- test/test_misc.f90 | 10 ++-- test/test_omega_prime.f90 | 85 ++++++++++++++--------------- test/test_torque.f90 | 21 ++++---- 16 files changed, 314 insertions(+), 306 deletions(-) diff --git a/src/attenuation_factor.f90 b/src/attenuation_factor.f90 index 1e258805..060ff1d7 100644 --- a/src/attenuation_factor.f90 +++ b/src/attenuation_factor.f90 @@ -42,7 +42,7 @@ subroutine init_attenuation_data() close (iunitarg) argmin = argvals(1) - argstepinv = 1.0_dp/(argvals(2) - argvals(1)) + argstepinv = 1.0_dp / (argvals(2) - argvals(1)) attenuation_data_initialized = .true. end subroutine init_attenuation_data diff --git a/src/collis_nbi.f90 b/src/collis_nbi.f90 index be22e44d..c715b73b 100644 --- a/src/collis_nbi.f90 +++ b/src/collis_nbi.f90 @@ -44,24 +44,24 @@ subroutine coleff(p, dpp, dhh, fpeff) fpeff = 0.0_dp do i = 1, nsorts - xbeta = p*velrat(i) + xbeta = p * velrat(i) call onseff(xbeta, d_p, dh, dpd) - dpp = dpp + d_p*efcolf(i) - dhh = dhh + dh*efcolf(i) - fpeff = fpeff + (dpd/plim - 2.0_dp*d_p*p*enrat(i))*efcolf(i) + dpp = dpp + d_p * efcolf(i) + dhh = dhh + dh * efcolf(i) + fpeff = fpeff + (dpd / plim - 2.0_dp * d_p * p * enrat(i)) * efcolf(i) end do - dhh = dhh/plim**2 + dhh = dhh / plim**2 return end subroutine coleff subroutine onseff(v, d_p, dh, dpd) - ! d_p - dimensionless dpp - ! dh - dhh*p^2 (p - dmls) - ! dpd - (1/p)(d/dp)p^2*d_p (p - dmls) + ! d_p - dimensionless dpp + ! dh - dhh*p^2 (p - dmls) + ! dpd - (1/p)(d/dp)p^2*d_p (p - dmls) implicit none @@ -74,21 +74,21 @@ subroutine onseff(v, d_p, dh, dpd) real(dp) :: v2, v3, ex, er v2 = v**2 - v3 = v2*v + v3 = v2 * v if (v < 0.01_dp) then - d_p = cons*(1.0_dp - 0.6_dp*v2) - dh = cons*(1.0_dp - 0.2_dp*v2) - dpd = 2.0_dp*cons*(1.0_dp - 1.2_dp*v2) + d_p = cons * (1.0_dp - 0.6_dp * v2) + dh = cons * (1.0_dp - 0.2_dp * v2) + dpd = 2.0_dp * cons * (1.0_dp - 1.2_dp * v2) elseif (v > 6.0_dp) then - d_p = 1.0_dp/v3 - dh = (1.0_dp - 0.5_dp/v2)/v - dpd = -1.0_dp/v3 + d_p = 1.0_dp / v3 + dh = (1.0_dp - 0.5_dp / v2) / v + dpd = -1.0_dp / v3 else - ex = exp(-v2)/sqp + ex = exp(-v2) / sqp er = erf(v) - d_p = er/v3 - 2.0_dp*ex/v2 - dh = er*(1.0_dp - 0.5_dp/v2)/v + ex/v2 - dpd = 4.0_dp*ex - d_p + d_p = er / v3 - 2.0_dp * ex / v2 + dh = er * (1.0_dp - 0.5_dp / v2) / v + ex / v2 + dpd = 4.0_dp * ex - d_p end if return @@ -138,33 +138,35 @@ subroutine loacol_nbi(amb, am1, am2, Zb, Z1, Z2, densi1, densi2, tempi1, tempi2, e = 4.8032e-10_dp ev = 1.6022e-12_dp - enrat(1) = ebeam/tempi1 - enrat(2) = ebeam/tempi2 - enrat(3) = ebeam/tempe + enrat(1) = ebeam / tempi1 + enrat(2) = ebeam / tempi2 + enrat(3) = ebeam / tempe - v0 = sqrt(2.0_dp*ebeam*ev/(amb*pmass)) - vti1 = sqrt(2.0_dp*tempi1*ev/(pmass*am1)) - vti2 = sqrt(2.0_dp*tempi2*ev/(pmass*am2)) - vte = sqrt(2.0_dp*tempe*ev/emass) + v0 = sqrt(2.0_dp * ebeam * ev / (amb * pmass)) + vti1 = sqrt(2.0_dp * tempi1 * ev / (pmass * am1)) + vti2 = sqrt(2.0_dp * tempi2 * ev / (pmass * am2)) + vte = sqrt(2.0_dp * tempe * ev / emass) - velrat(1) = v0/vti1 - velrat(2) = v0/vti2 - velrat(3) = v0/vte + velrat(1) = v0 / vti1 + velrat(2) = v0 / vti2 + velrat(3) = v0 / vte - dense = densi1*Z1 + densi2*Z2 + dense = densi1 * Z1 + densi2 * Z2 alami1 = 23.0_dp - log(max(epsilon(1.0_dp), & - sqrt(densi1*Z1**2/tempi1)*Zb*Z1*(amb + am1)/(amb*tempi1 + am1*ebeam))) + sqrt(densi1 * Z1**2 / tempi1) * Zb * Z1 * (amb + am1) / (amb * & + tempi1 + am1 * ebeam))) alami2 = 23.0_dp - log(max(epsilon(1.0_dp), & - sqrt(densi2*Z2**2/tempi2)*Zb*Z2*(amb + am2)/(amb*tempi2 + am2*ebeam))) - alame = 24.0_dp - log(sqrt(dense)/tempe) - frecol_base = 2.0_dp*pi*dense*e**4*Zb**2/((amb*pmass)**2*v0**3) ! usual - frecol_base = frecol_base/v0 ! normalized + sqrt(densi2 * Z2**2 / tempi2) * Zb * Z2 * (amb + am2) / (amb * & + tempi2 + am2 * ebeam))) + alame = 24.0_dp - log(sqrt(dense) / tempe) + frecol_base = 2.0_dp * pi * dense * e**4 * Zb**2 / ((amb * pmass)**2 * v0**3) ! usual + frecol_base = frecol_base / v0 ! normalized - efcolf(1) = frecol_base*Z1**2*alami1*densi1/dense - efcolf(2) = frecol_base*Z2**2*alami2*densi2/dense - efcolf(3) = frecol_base*alame + efcolf(1) = frecol_base * Z1**2 * alami1 * densi1 / dense + efcolf(2) = frecol_base * Z2**2 * alami2 * densi2 / dense + efcolf(3) = frecol_base * alame - efcolf = efcolf*velrat + efcolf = efcolf * velrat end subroutine loacol_nbi end module collis_alp diff --git a/src/do_magfie_standalone.f90 b/src/do_magfie_standalone.f90 index 4aab6923..fcb3f714 100644 --- a/src/do_magfie_standalone.f90 +++ b/src/do_magfie_standalone.f90 @@ -6,7 +6,7 @@ module do_magfie_mod implicit none real(dp), private :: s_prev = -1.0_dp - real(dp), private, allocatable :: spl_val_c(:,:), spl_val_s(:,:) + real(dp), private, allocatable :: spl_val_c(:, :), spl_val_s(:, :) ! Work arrays previously automatic on stack (size=nmode) real(dp), private, allocatable :: B0mnc(:), dB0dsmnc(:), B0mns(:), dB0dsmns(:) real(dp), private, allocatable :: costerm(:), sinterm(:) @@ -14,9 +14,9 @@ module do_magfie_mod real(dp), parameter :: sign_theta = -1.0_dp ! negative for left-handed real(dp) :: s = 0.0_dp, psi_pr = 0.0_dp, Bthcov = 0.0_dp, Bphcov = 0.0_dp, & - dBthcovds = 0.0_dp, dBphcovds = 0.0_dp, & - q = 0.0_dp, dqds = 0.0_dp, iota = 0.0_dp, R0 = 0.0_dp, a = 0.0_dp, & - eps = 0.0_dp, B0h = 0.0_dp, B00 = 0.0_dp + dBthcovds = 0.0_dp, dBphcovds = 0.0_dp, & + q = 0.0_dp, dqds = 0.0_dp, iota = 0.0_dp, R0 = 0.0_dp, a = 0.0_dp, & + eps = 0.0_dp, B0h = 0.0_dp, B00 = 0.0_dp real(dp) :: bfac = 1.0_dp ! B0h is the 0th theta harmonic of bmod on current flux surface ! and B00 the 0th theta harmonic of bmod on the innermost flux surface @@ -28,7 +28,7 @@ module do_magfie_mod ! Work arrays for booz_to_cyl (size=nmode) real(dp), private, allocatable :: rmnc(:), rmns(:), zmnc(:), zmns(:) - real(dp), parameter :: ItoB = 2.0e-1_dp*sign_theta ! Covarient B (cgs) from I (SI) + real(dp), parameter :: ItoB = 2.0e-1_dp * sign_theta ! Covarient B (cgs) from I (SI) ! Bcov=mu0/2pi*I,mu0->4pi/c,I->10^(-1)*c*I integer :: ncol1 = 0, ncol2 = 0 ! number of columns in input file @@ -104,7 +104,7 @@ subroutine read_boozer_file(path) end do ! Set B00 from first mode - B00 = 1.0e4_dp*modes0(1, 1, 6)*bfac + B00 = 1.0e4_dp * modes0(1, 1, 6) * bfac end subroutine read_boozer_file subroutine init_magfie_at_s() @@ -187,56 +187,56 @@ subroutine do_magfie(x, bmod, sqrtg, bder, hcovar, hctrvr, hcurl) x1 = min(params0(nflux, 1), x1) spl_val = spline_val_0(spl_coeff1(:, :, 3), x1) - Bthcov = ItoB*spl_val(1)*bfac - dBthcovds = ItoB*spl_val(2)*bfac + Bthcov = ItoB * spl_val(1) * bfac + dBthcovds = ItoB * spl_val(2) * bfac spl_val = spline_val_0(spl_coeff1(:, :, 2), x1) - Bphcov = ItoB*spl_val(1)*bfac - dBphcovds = ItoB*spl_val(2)*bfac + Bphcov = ItoB * spl_val(1) * bfac + dBphcovds = ItoB * spl_val(2) * bfac spl_val = spline_val_0(spl_coeff1(:, :, 1), x1) iota = spl_val(1) - q = 1/iota - dqds = -spl_val(2)/iota**2 + q = 1 / iota + dqds = -spl_val(2) / iota**2 call fast_sin_cos(modes0(1, :, 1), x(3), sinterm, costerm) ! calculate B-field from modes if (inp_swi == 8) then call cached_spline(x1, s_prev, spl_coeff2(:, :, 4, :), spl_val_c) - B0mnc(:) = 1.0e4_dp*spl_val_c(1, :)*bfac - dB0dsmnc(:) = 1.0e4_dp*spl_val_c(2, :)*bfac + B0mnc(:) = 1.0e4_dp * spl_val_c(1, :) * bfac + dB0dsmnc(:) = 1.0e4_dp * spl_val_c(2, :) * bfac B0h = B0mnc(1) - bmod = sum(B0mnc*costerm) - bder(1) = sum(dB0dsmnc*costerm)/bmod + bmod = sum(B0mnc * costerm) + bder(1) = sum(dB0dsmnc * costerm) / bmod bder(2) = 0.0_dp - bder(3) = sum(-modes0(1, :, 1)*B0mnc*sinterm)/bmod + bder(3) = sum(-modes0(1, :, 1) * B0mnc * sinterm) / bmod else if (inp_swi == 9) then call cached_spline(x1, s_prev, spl_coeff2(:, :, 7, :), spl_val_c) - B0mnc(:) = 1.0e4_dp*spl_val_c(1, :)*bfac - dB0dsmnc(:) = 1.0e4_dp*spl_val_c(2, :)*bfac + B0mnc(:) = 1.0e4_dp * spl_val_c(1, :) * bfac + dB0dsmnc(:) = 1.0e4_dp * spl_val_c(2, :) * bfac call cached_spline(x1, s_prev, spl_coeff2(:, :, 8, :), spl_val_s) - B0mns(:) = 1.0e4_dp*spl_val_s(1, :)*bfac - dB0dsmns(:) = 1.0e4_dp*spl_val_s(2, :)*bfac + B0mns(:) = 1.0e4_dp * spl_val_s(1, :) * bfac + dB0dsmns(:) = 1.0e4_dp * spl_val_s(2, :) * bfac B0h = B0mnc(1) - bmod = sum(B0mnc*costerm + B0mns*sinterm) - bder(1) = sum(dB0dsmnc*costerm + dB0dsmns*sinterm)/bmod + bmod = sum(B0mnc * costerm + B0mns * sinterm) + bder(1) = sum(dB0dsmnc * costerm + dB0dsmns * sinterm) / bmod bder(2) = 0.0_dp - bder(3) = sum(-modes0(1, :, 1)*B0mnc*sinterm & - + modes0(1, :, 1)*B0mns*costerm)/bmod + bder(3) = sum(-modes0(1, :, 1) * B0mnc * sinterm & + + modes0(1, :, 1) * B0mns * costerm) / bmod end if - sqgbmod2 = sign_theta*psi_pr*(Bphcov + iota*Bthcov) - sqgbmod = sqgbmod2/bmod - sqrtg = sqgbmod/bmod + sqgbmod2 = sign_theta * psi_pr * (Bphcov + iota * Bthcov) + sqgbmod = sqgbmod2 / bmod + sqrtg = sqgbmod / bmod hcovar(1) = 0.0_dp ! TODO - hcovar(2) = Bphcov/bmod - hcovar(3) = Bthcov/bmod + hcovar(2) = Bphcov / bmod + hcovar(3) = Bthcov / bmod hctrvr(1) = 0.0_dp - hctrvr(2) = sign_theta*psi_pr/sqgbmod - hctrvr(3) = sign_theta*iota*psi_pr/sqgbmod + hctrvr(2) = sign_theta * psi_pr / sqgbmod + hctrvr(3) = sign_theta * iota * psi_pr / sqgbmod hcurl(1) = 0.0_dp ! TODO hcurl(3) = 0.0_dp ! TODO @@ -255,12 +255,12 @@ subroutine boozer_read(filename) open (unit=18, file=filename, action='read', status='old') read (18, '(////)') read (18, *) m0b, n0b, nflux, nfp, flux, a, R0 - a = 100*a ! m -> cm - R0 = 100*R0 ! m -> cm + a = 100 * a ! m -> cm + R0 = 100 * R0 ! m -> cm - psi_pr = 1.0e8_dp*flux/(2*pi)*bfac ! T -> Gauss, m -> cm + psi_pr = 1.0e8_dp * flux / (2 * pi) * bfac ! T -> Gauss, m -> cm - nmode = (m0b + 1)*(n0b + 1) + nmode = (m0b + 1) * (n0b + 1) ! Allocate params and modes (deallocate first if size changed) if (allocated(params0)) then @@ -299,18 +299,18 @@ subroutine booz_to_cyl(x, r) do j = 1, nmode spl_val = spline_val_0(spl_coeff2(:, :, 1, j), x1) - rmnc(j) = 1.0e2_dp*spl_val(1) + rmnc(j) = 1.0e2_dp * spl_val(1) spl_val = spline_val_0(spl_coeff2(:, :, 2, j), x1) - rmns(j) = 1.0e2_dp*spl_val(1) + rmns(j) = 1.0e2_dp * spl_val(1) spl_val = spline_val_0(spl_coeff2(:, :, 3, j), x1) - zmnc(j) = 1.0e2_dp*spl_val(1) + zmnc(j) = 1.0e2_dp * spl_val(1) spl_val = spline_val_0(spl_coeff2(:, :, 4, j), x1) - zmns(j) = 1.0e2_dp*spl_val(1) + zmns(j) = 1.0e2_dp * spl_val(1) end do - r(1) = sum(rmnc*cos(modes0(1, :, 1)*x(3)) + rmns*sin(modes0(1, :, 1)*x(3))) + r(1) = sum(rmnc * cos(modes0(1, :, 1) * x(3)) + rmns * sin(modes0(1, :, 1) * x(3))) r(2) = 0.0_dp ! TODO: phi - r(3) = sum(zmnc*cos(modes0(1, :, 1)*x(3)) + zmns*sin(modes0(1, :, 1)*x(3))) + r(3) = sum(zmnc * cos(modes0(1, :, 1) * x(3)) + zmns * sin(modes0(1, :, 1) * x(3))) end subroutine booz_to_cyl @@ -347,7 +347,7 @@ module do_magfie_pert_mod implicit none real(dp), private :: s_prev = -1.0_dp - real(dp), private, allocatable :: spl_val_c(:,:), spl_val_s(:,:) + real(dp), private, allocatable :: spl_val_c(:, :), spl_val_s(:, :) real(dp), allocatable, protected :: params(:, :), modes(:, :, :) integer, protected :: mb, nb, nflux, nfp, nmode @@ -357,9 +357,9 @@ module do_magfie_pert_mod ! Work arrays (size=nmode) real(dp), private, allocatable :: Bmnc(:), Bmns(:) - integer :: ncol1, ncol2 ! number of columns in input file - real(dp) :: mph ! toroidal perturbation mode (threadprivate) - real(dp) :: mph_shared = 0.0_dp ! shared copy for namelist input (when pertfile=.false.) + integer :: ncol1, ncol2 ! number of columns in input file + real(dp) :: mph ! toroidal perturbation mode (threadprivate) + real(dp) :: mph_shared = 0.0_dp ! shared copy for namelist input (when pertfile=.false.) ! Initialization flag for threadprivate allocatable arrays logical, save :: magfie_pert_arrays_initialized = .false. @@ -496,13 +496,13 @@ subroutine do_magfie_pert_amp(x, bamp) ! calculate B-field from modes if (inp_swi == 8) then call cached_spline(x1, s_prev, spl_coeff2(:, :, 4, :), spl_val_c) - Bmnc(:) = 1.0e4_dp*spl_val_c(1, :)*bfac - bamp = sum(Bmnc*cos(modes(1, :, 1)*x(3))) + Bmnc(:) = 1.0e4_dp * spl_val_c(1, :) * bfac + bamp = sum(Bmnc * cos(modes(1, :, 1) * x(3))) else if (inp_swi == 9) then call cached_spline(x1, s_prev, spl_coeff2(:, :, 7, :), spl_val_c) call cached_spline(x1, s_prev, spl_coeff2(:, :, 8, :), spl_val_s) - Bmnc(:) = 1.0e4_dp*spl_val_c(1, :)*bfac - Bmns(:) = 1.0e4_dp*spl_val_s(1, :)*bfac + Bmnc(:) = 1.0e4_dp * spl_val_c(1, :) * bfac + Bmns(:) = 1.0e4_dp * spl_val_s(1, :) * bfac bamp = fast_fourier_sum(Bmnc, Bmns, modes(1, :, 1), x(3)) end if diff --git a/src/freq.f90 b/src/freq.f90 index e4a9c0f0..874a3a7e 100644 --- a/src/freq.f90 +++ b/src/freq.f90 @@ -85,9 +85,9 @@ subroutine init_canon_freq_trapped_spline OmtB1 = 0.0_dp v = vth - etamin = (1.0_dp + epst)*etatp - etamax = etatp + (etadt - etatp)*(1.0_dp - epsst_spl) - ! Allocate coefficient arrays for trapped region splines (safe for undefined allocation status) + etamin = (1.0_dp + epst) * etatp + etamax = etatp + (etadt - etatp) * (1.0_dp - epsst_spl) + ! Allocate coefficient arrays for trapped region splines (safe for undefined allocation status) if (.not. freq_trapped_initialized) then if (allocated(Omth_spl_coeff)) deallocate(Omth_spl_coeff) if (allocated(OmtB_spl_coeff)) deallocate(OmtB_spl_coeff) @@ -109,10 +109,10 @@ subroutine init_canon_freq_trapped_spline ! logspace for eta b = log(epst_spl) - aa = 1.0_dp/(netaspl - 1.0_dp)*(log(etamax/etamin - 1.0_dp) - b) + aa = 1.0_dp / (netaspl - 1.0_dp) * (log(etamax / etamin - 1.0_dp) - b) do k = netaspl - 1, 0, -1 - eta = etamin*(1.0_dp + exp(aa*k + b)) + eta = etamin * (1.0_dp + exp(aa * k + b)) etarange(k + 1) = eta if (get_log_level() >= LOG_TRACE) then write(*,'(A,I4,A,ES12.5)') '[TRACE] init_canon_freq_trapped_spline k=', k, ' eta=', eta @@ -202,11 +202,11 @@ subroutine init_canon_freq_passing_spline write(*,'(A,1X,I0,2X,A,1X,I0,2X,A,1X,L1)') ' mph =', int(mph), 'sign_vpar =', int(sign_vpar), 'nonlin =', nonlin end if - b = log((etamax - etamin)/etamax) - aa = 1.0_dp/(netaspl_pass - 1.0_dp)*(log(epsp_spl) - b) + b = log((etamax - etamin) / etamax) + aa = 1.0_dp / (netaspl_pass - 1.0_dp) * (log(epsp_spl) - b) do k = netaspl_pass - 1, 0, -1 - eta = etamax*(1.0_dp - exp(aa*k + b)) + eta = etamax * (1.0_dp - exp(aa * k + b)) etarange(k + 1) = eta if (get_log_level() >= LOG_TRACE) then write(*,'(A,I4,A,ES12.5)') '[TRACE] init_canon_freq_passing_spline k=', k, ' eta=', eta @@ -255,27 +255,29 @@ subroutine Om_tB(v, eta, OmtB, dOmtBdv, dOmtBdeta) real(dp) :: splineval(3) real(dp) :: Omth, dOmthdv, dOmthdeta if (eta > etatp) then - if (eta > etatp*(1 + epst_spl)) then + if (eta > etatp * (1 + epst_spl)) then splineval = spline_val_0(OmtB_spl_coeff, eta) else ! extrapolation call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) - splineval(1) = sign_vpar*(k_OmtB_t*log(eta - etatp) + d_OmtB_t)*Omth/v - splineval(2) = sign_vpar*(Omth/v*k_OmtB_t/(eta - etatp) + & - dOmthdeta/v*(k_OmtB_t*log(eta - etatp) + d_OmtB_t)) + splineval(1) = sign_vpar * (k_OmtB_t * log(eta - etatp) + d_OmtB_t) * Omth / v + splineval(2) = sign_vpar * (Omth / v * k_OmtB_t / (eta - etatp) + & + dOmthdeta / v * (k_OmtB_t * log(eta - etatp) & + + d_OmtB_t)) end if else - if (eta < etatp*(1 - epsp_spl)) then + if (eta < etatp * (1 - epsp_spl)) then splineval = spline_val_0(OmtB_pass_spl_coeff, eta) else ! extrapolation call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) - splineval(1) = sign_vpar*(k_OmtB_p*log(etatp - eta) + d_OmtB_p)*Omth/v - splineval(2) = sign_vpar*(Omth/v*k_OmtB_p/(eta - etatp) + & - dOmthdeta/v*(k_OmtB_p*log(etatp - eta) + d_OmtB_p)) + splineval(1) = sign_vpar * (k_OmtB_p * log(etatp - eta) + d_OmtB_p) * Omth / v + splineval(2) = sign_vpar * (Omth / v * k_OmtB_p / (eta - etatp) + & + dOmthdeta / v * (k_OmtB_p * log(etatp - eta) & + + d_OmtB_p)) end if end if - OmtB = splineval(1)*v**2 - dOmtBdv = 2.0_dp*splineval(1)*v - dOmtBdeta = splineval(2)*v**2 + OmtB = splineval(1) * v**2 + dOmtBdv = 2.0_dp * splineval(1) * v + dOmtBdeta = splineval(2) * v**2 end subroutine Om_tB subroutine Om_ph(v, eta, Omph, dOmphdv, dOmphdeta) @@ -320,21 +322,21 @@ subroutine Om_th(v, eta, Omth, dOmthdv, dOmthdeta) if (eta > etatp) then if (eta > etatp*(1 + epst_spl)) then splineval = spline_val_0(Omth_spl_coeff, eta) - else ! extrapolation - splineval(1) = 2.0_dp*pi/(k_taub_t*log(eta - etatp) + d_taub_t) - splineval(2) = -splineval(1)**2/(2.0_dp*pi)*k_taub_t/(eta - etatp) + else ! extrapolation + splineval(1) = 2.0_dp * pi / (k_taub_t * log(eta - etatp) + d_taub_t) + splineval(2) = -splineval(1)**2 / (2.0_dp * pi) * k_taub_t / (eta - etatp) end if else - if (eta < etatp*(1 - epsp_spl)) then + if (eta < etatp * (1 - epsp_spl)) then splineval = spline_val_0(Omth_pass_spl_coeff, eta) - else ! extrapolation - splineval(1) = 2.0_dp*pi/(k_taub_p*log(etatp - eta) + d_taub_p) - splineval(2) = -splineval(1)**2/(2.0_dp*pi)*k_taub_p/(eta - etatp) + else ! extrapolation + splineval(1) = 2.0_dp * pi / (k_taub_p * log(etatp - eta) + d_taub_p) + splineval(2) = -splineval(1)**2 / (2.0_dp * pi) * k_taub_p / (eta - etatp) end if end if - Omth = sign_vpar*splineval(1)*v - dOmthdv = sign_vpar*splineval(1) - dOmthdeta = sign_vpar*splineval(2)*v + Omth = sign_vpar * splineval(1) * v + dOmthdv = sign_vpar * splineval(1) + dOmthdeta = sign_vpar * splineval(2) * v end subroutine Om_th subroutine d_Om_ds(v, eta, taub_estimate, dOmthds, dOmphds) @@ -349,40 +351,41 @@ subroutine d_Om_ds(v, eta, taub_estimate, dOmthds, dOmphds) s0 = s ds = 2.0e-8_dp - s = s0 - ds/2.0_dp + s = s0 - ds / 2.0_dp taub_est = bounce_time(v, eta, taub_estimate) taub = taub_est call bounce_fast(v, eta, taub, bounceavg, timestep) - Omth = sign_vpar_htheta*2.0_dp*pi/taub + Omth = sign_vpar_htheta * 2.0_dp * pi / taub if (magdrift) then if (eta > etatp) then - Omph_noE = bounceavg(3)*v**2 + Omph_noE = bounceavg(3) * v**2 else - Omph_noE = bounceavg(3)*v**2 + Omth/iota + Omph_noE = bounceavg(3) * v**2 + Omth / iota end if else if (eta > etatp) then Omph_noE = 0.0_dp else - Omph_noE = Omth/iota + Omph_noE = Omth / iota end if end if - s = s0 + ds/2.0_dp + s = s0 + ds / 2.0_dp taub_est = bounce_time(v, eta, taub_estimate) taub = taub_est call bounce_fast(v, eta, taub, bounceavg, timestep) - dOmthds = sign_vpar_htheta*(2.0_dp*pi/taub - sign_vpar_htheta*Omth)/ds + dOmthds = sign_vpar_htheta * (2.0_dp * pi / taub - sign_vpar_htheta * Omth) / ds if (magdrift) then if (eta > etatp) then - dOmphds = dOm_tEds + (bounceavg(3)*v**2 - Omph_noE)/ds + dOmphds = dOm_tEds + (bounceavg(3) * v**2 - Omph_noE) / ds else - dOmphds = dOm_tEds + (bounceavg(3)*v**2 + (2.0_dp*pi/taub)/iota - Omph_noE)/ds + dOmphds = dOm_tEds + (bounceavg(3) * v**2 + (2.0_dp * pi / taub) / iota - & + Omph_noE) / ds end if else if (eta > etatp) then dOmphds = dOm_tEds else - dOmphds = dOm_tEds + ((2.0_dp*pi/taub)/iota - Omph_noE)/ds + dOmphds = dOm_tEds + ((2.0_dp * pi / taub) / iota - Omph_noE) / ds end if end if diff --git a/src/magfie.f90 b/src/magfie.f90 index 87f92c9a..27193dae 100644 --- a/src/magfie.f90 +++ b/src/magfie.f90 @@ -54,12 +54,12 @@ subroutine init_flux_surface_average(s) if (bmod > Bmax) Bmax = bmod end do - dVds = 2.0_dp*pi*dVds - B0 = B0/(2.0_dp*pi) - eps = eps/(B0*pi) + dVds = 2.0_dp * pi * dVds + B0 = B0 / (2.0_dp * pi) + eps = eps / (B0 * pi) - etatp = 1.0_dp/Bmax - etadt = 1.0_dp/Bmin + etatp = 1.0_dp / Bmax + etadt = 1.0_dp / Bmin write(buffer, "(A,ES12.5)") " eps calc: ", eps call log_result(buffer) diff --git a/src/neort.f90 b/src/neort.f90 index 3531e065..a31276ea 100644 --- a/src/neort.f90 +++ b/src/neort.f90 @@ -65,8 +65,8 @@ subroutine check_magfie(data) real(dp) :: dpp, dhh, fpeff complex(dp) :: bn - D_plateau = pi*vth**3/(16.0_dp*R0*iota*(qi*B0/(mi*c))**2) - Drp = 4*mph*q/(eps**2*sqrt(pi)) ! actually this is Drp/Dp + D_plateau = pi * vth**3 / (16.0_dp * R0 * iota * (qi * B0 / (mi * c))**2) + Drp = 4 * mph * q / (eps**2 * sqrt(pi)) ! actually this is Drp/Dp thmin = -pi thmax = pi @@ -89,11 +89,11 @@ subroutine check_magfie(data) data%params%dVds = dVds data%params%M_t = M_t data%params%Om_tE = Om_tE - data%params%Om_tBref = c*mi*vth**2/(2.0_dp*qi*sign_theta*psi_pr) + data%params%Om_tBref = c * mi * vth**2 / (2.0_dp * qi * sign_theta * psi_pr) data%params%vth = vth - data%params%T_in_eV = mi/2.0_dp*vth**2/eV - data%params%m0 = 1.0_dp*m0 - data%params%n0 = 1.0_dp*mph + data%params%T_in_eV = mi / 2.0_dp * vth**2 / eV + data%params%m0 = 1.0_dp * m0 + data%params%n0 = 1.0_dp * mph data%params%Dp = D_plateau data%params%Drp_over_Dp = Drp data%params%etatp = etatp @@ -227,8 +227,8 @@ subroutine compute_transport_harmonic(j, Dco, Dctr, Dt, Tco, Tctr, Tt, harmonic) Tresctr = 0.0_dp Trest = 0.0_dp - vminp = 1.0e-6_dp*vth - vmaxp = 3.0_dp*vth + vminp = 1.0e-6_dp * vth + vmaxp = 3.0_dp * vth vmint = vminp vmaxt = vmaxp diff --git a/src/nonlin.f90 b/src/nonlin.f90 index 7a20bf0d..78b15c72 100644 --- a/src/nonlin.f90 +++ b/src/nonlin.f90 @@ -26,20 +26,20 @@ function nonlinear_attenuation(ux, eta, bounceavg, Omth, dOmthdv, dOmthdeta, Hmn return end if - v = ux*vth + v = ux * vth if (nonlin) then call Om_ph(v, eta, Omph, dOmphdv, dOmphdeta) - call d_Om_ds(v, eta, 2.0_dp*pi/Omth, dOmthds, dOmphds) - dOmdv = mth*dOmthdv + mph*dOmphdv - dOmdeta = mth*dOmthdeta + mph*dOmphdeta - dOmdpph = -(qi/c*iota*sign_theta*psi_pr)**(-1)*(mth*dOmthds + mph*dOmphds) + call d_Om_ds(v, eta, 2.0_dp * pi / Omth, dOmthds, dOmphds) + dOmdv = mth * dOmthdv + mph * dOmphdv + dOmdeta = mth * dOmthdeta + mph * dOmphdeta + dOmdpph = -(qi / c * iota * sign_theta * psi_pr)**(-1) * (mth * dOmthds + mph * dOmphds) Ompr = omega_prime(ux, eta, bounceavg, Omth, dOmdv, dOmdeta, dOmdpph) call coleff(ux, dpp, dhh, fpeff) - dhh = vth*dhh - dpp = vth**3*dpp - dres = dpp*(dOmdv/Ompr)**2 + dhh*eta*(bounceavg(5) - eta)*(dOmdeta/Ompr)**2 - dnorm = dres*sqrt(abs(Ompr))/sqrt(abs(Hmn2))**(3.0_dp/2.0_dp) + dhh = vth * dhh + dpp = vth**3 * dpp + dres = dpp * (dOmdv / Ompr)**2 + dhh * eta * (bounceavg(5) - eta) * (dOmdeta / Ompr)**2 + dnorm = dres * sqrt(abs(Ompr)) / sqrt(abs(Hmn2))**(3.0_dp / 2.0_dp) call attenuation_factor(dnorm, nonlinear_attenuation) end if @@ -52,17 +52,17 @@ pure function omega_prime(ux, eta, bounceavg, Omth, dOmdv, dOmdeta, dOmdpph) real(dp) :: ma, mb, mc, md, me, mf, dvdJ, detadJ - ma = mi*(ux*vth)*mi*c/qi*eta - mb = mi*(ux*vth)**2/2*mi*c/qi - mc = mi/(2.0_dp*Omth)*(ux*vth)*(1.0_dp - eta*bounceavg(6)) - md = mi*(ux*vth)**2/2.0_dp*Omth - me = -mth/mph - mf = 1.0_dp/mph + ma = mi * (ux * vth) * mi * c / qi * eta + mb = mi * (ux * vth)**2 / 2 * mi * c / qi + mc = mi / (2.0_dp * Omth) * (ux * vth) * (1.0_dp - eta * bounceavg(6)) + md = mi * (ux * vth)**2 / 2.0_dp * Omth + me = -mth / mph + mf = 1.0_dp / mph - dvdJ = mb*me/(ma*md*mf - mb*mc*mf) - detadJ = ma*me/(mb*mc*mf - ma*md*mf) + dvdJ = mb * me / (ma * md * mf - mb * mc * mf) + detadJ = ma * me / (mb * mc * mf - ma * md * mf) - omega_prime = dOmdv*dvdJ + dOmdeta*detadJ + mph*dOmdpph + omega_prime = dOmdv * dvdJ + dOmdeta * detadJ + mph * dOmdpph end function omega_prime end module neort_nonlin diff --git a/src/orbit.f90 b/src/orbit.f90 index 9f2cf66b..713bfa99 100644 --- a/src/orbit.f90 +++ b/src/orbit.f90 @@ -87,12 +87,12 @@ subroutine bounce(v, eta, taub, bounceavg, taub_estimate) call evaluate_bfield_local(bmod, htheta) sign_vpar_htheta = sign(1.0_dp, htheta)*sign_vpar y0 = 1.0e-15_dp - y0(1) = th0 ! poloidal angle theta - y0(2) = sign_vpar_htheta*vpar(v, eta, bmod) ! parallel velocity vpar - y0(3) = 0.0_dp ! toroidal velocity v_ph for drift frequency Om_ph - y0(4) = 0.0_dp ! perturbed Hamiltonian real part - y0(5) = 0.0_dp ! perturbed Hamiltonian imaginary part - y0(6) = 0.0_dp ! 1/abs(B) + y0(1) = th0 ! poloidal angle theta + y0(2) = sign_vpar_htheta * vpar(v, eta, bmod) ! parallel velocity vpar + y0(3) = 0.0_dp ! toroidal velocity v_ph for drift frequency Om_ph + y0(4) = 0.0_dp ! perturbed Hamiltonian real part + y0(5) = 0.0_dp ! perturbed Hamiltonian imaginary part + y0(6) = 0.0_dp ! 1/abs(B) ! y0(7) = 0.0_dp ! abs(B) ! If bounce time estimate exists (elliptic integrals), @@ -100,15 +100,15 @@ subroutine bounce(v, eta, taub, bounceavg, taub_estimate) if (present(taub_estimate)) then taub = taub_estimate else - taub = 2.0*pi/abs(vperp(v, eta, bmod)*iota/R0*sqrt(eps/2.0_dp)) + taub = 2.0 * pi / abs(vperp(v, eta, bmod) * iota / R0 * sqrt(eps / 2.0_dp)) end if ! Look for exactly one orbit turn via root-finding. ! Start by looking for 5 points per turn. - findroot_res = bounce_integral(v, eta, nvar, y0, taub/5.0_dp, timestep) + findroot_res = bounce_integral(v, eta, nvar, y0, taub / 5.0_dp, timestep) taub = findroot_res(1) - bounceavg = findroot_res(2:)/taub + bounceavg = findroot_res(2:) / taub end subroutine bounce subroutine evaluate_bfield_local(bmod, htheta) @@ -126,7 +126,7 @@ pure function vpar(v, eta, bmod) ! parallel velocity real(dp) :: vpar real(dp), intent(in) :: v, eta, bmod - vpar = v*sqrt(1.0_dp - eta*bmod) + vpar = v * sqrt(1.0_dp - eta * bmod) if (isnan(vpar)) then vpar = 0.0_dp end if @@ -162,10 +162,10 @@ subroutine bounce_fast(v, eta, taub, bounceavg, ts, istate_out) t2 = taub call evaluate_bfield_local(bmod, htheta) - sign_vpar_htheta = sign(1.0_dp, htheta)*sign_vpar + sign_vpar_htheta = sign(1.0_dp, htheta) * sign_vpar y = 1.0e-15_dp y(1) = th0 - y(2) = sign_vpar_htheta*vpar(v, eta, bmod) + y(2) = sign_vpar_htheta * vpar(v, eta, bmod) y(3:6) = 0.0_dp neq = nvar @@ -205,20 +205,20 @@ function bounce_time(v, eta, taub_estimate) result(taub) real(dp) :: taub integer, parameter :: neq = 2 - real(dp) :: y0(neq), roots(neq+1) + real(dp) :: y0(neq), roots(neq + 1) real(dp) :: bmod, htheta call trace('bounce_time') call evaluate_bfield_local(bmod, htheta) - sign_vpar_htheta = sign(1.0_dp, htheta)*sign_vpar + sign_vpar_htheta = sign(1.0_dp, htheta) * sign_vpar - y0(1) = th0 ! poloidal angle theta - y0(2) = sign_vpar_htheta*vpar(v, eta, bmod) ! parallel velocity vpar + y0(1) = th0 ! poloidal angle theta + y0(2) = sign_vpar_htheta * vpar(v, eta, bmod) ! parallel velocity vpar if (present(taub_estimate)) then taub = taub_estimate else - taub = 2.0*pi/abs(vperp(v, eta, bmod)*iota/R0*sqrt(eps/2.0_dp)) + taub = 2.0 * pi / abs(vperp(v, eta, bmod) * iota / R0 * sqrt(eps / 2.0_dp)) end if roots = bounce_integral(v, eta, neq, y0, taub, timestep_poloidal_motion) @@ -247,8 +247,8 @@ pure subroutine poloidal_velocity(v, eta, bmod, hthctr, hderth, v_par, ydot) real(dp), intent(in) :: v_par real(dp), intent(out) :: ydot(2) - ydot(1) = v_par*hthctr ! theta - ydot(2) = -v**2*eta/2.0_dp*hthctr*hderth*bmod ! v_par + ydot(1) = v_par * hthctr ! theta + ydot(2) = -v**2 * eta / 2.0_dp * hthctr * hderth * bmod ! v_par end subroutine poloidal_velocity function bounce_integral(v, eta, neq, y0, dt, ts) @@ -353,8 +353,8 @@ subroutine bounceroots(NEQ, T, Y, NG, GOUT) real(dp), intent(out) :: GOUT(ng) associate (dummy => T) end associate - GOUT(1) = sign_vpar_htheta*(Y(1) - th0) ! trapped orbit return to starting point - GOUT(2) = sign_vpar_htheta*(2.0_dp*pi - (Y(1) - th0)) ! passing orbit return + GOUT(1) = sign_vpar_htheta * (Y(1) - th0) ! trapped orbit return to starting point + GOUT(2) = sign_vpar_htheta * (2.0_dp * pi - (Y(1) - th0)) ! passing orbit return return end subroutine bounceroots @@ -375,26 +375,25 @@ subroutine timestep(v, eta, neq, t, y, ydot) real(dp) :: Om_tB_v real(dp) :: shearterm - x(1) = s x(2) = 0.0_dp x(3) = y(1) call do_magfie(x, bmod, sqrtg, hder, hcovar, hctrvr, hcurl) - shearterm = Bphcov*dqds + shearterm = Bphcov * dqds if (noshear) then shearterm = 0 end if - Om_tB_v = mi*c*q/(2.0_dp*qi*sign_theta*psi_pr*bmod)*( & ! Om_tB/v**2 - -(2.0_dp - eta*bmod)*bmod*hder(1) & - + 2.0_dp*(1.0_dp - eta*bmod)*hctrvr(3)* & - (dBthcovds + q*dBphcovds + shearterm)) + Om_tB_v = mi * c * q / (2.0_dp * qi * sign_theta * psi_pr * bmod) * ( & ! Om_tB/v**2 + -(2.0_dp - eta * bmod) * bmod * hder(1) & + + 2.0_dp * (1.0_dp - eta * bmod) * hctrvr(3) * & + (dBthcovds + q * dBphcovds + shearterm)) - ydot(1) = y(2)*hctrvr(3) ! theta - ydot(2) = -0.5_dp*v**2*eta*hctrvr(3)*hder(3)*bmod ! v_par + ydot(1) = y(2) * hctrvr(3) ! theta + ydot(2) = -0.5_dp * v**2 * eta * hctrvr(3) * hder(3) * bmod ! v_par ydot(3) = Om_tB_v ! for bounce average of Om_tB/v**2 - ydot(4:) = 0.0_dp ! remaining integrands not computed here + ydot(4:) = 0.0_dp ! remaining integrands not computed here end subroutine timestep end module neort_orbit diff --git a/src/profiles.f90 b/src/profiles.f90 index fd046cfc..65b2a6ba 100644 --- a/src/profiles.f90 +++ b/src/profiles.f90 @@ -23,10 +23,10 @@ module neort_profiles ! Plasma profile spline data integer :: nplasma_global = 0 real(dp) :: am1_global = 0.0_dp, am2_global = 0.0_dp, Z1_global = 0.0_dp, Z2_global = 0.0_dp - real(dp), allocatable :: plasma_spl_coeff(:,:,:) + real(dp), allocatable :: plasma_spl_coeff(:, :, :) ! Rotation profile spline data - real(dp), allocatable :: Mt_spl_coeff(:,:) + real(dp), allocatable :: Mt_spl_coeff(:, :) ! Flux-surface dependent quantities (interpolated per-thread at each s) !$omp threadprivate (vth, dvthds, M_t, dM_tds, Om_tE, dOm_tEds) !$omp threadprivate (ni1, ni2, Ti1, Ti2, Te, dni1ds, dni2ds, dTi1ds, dTi2ds, dTeds) @@ -96,16 +96,16 @@ subroutine init_plasma_at_s() dTeds = spl_val(2) ! Compute derived quantities - qi = Z1_global*qe - mi = am1_global*mu - vth = sqrt(2.0_dp*Ti1*ev/mi) - dvthds = 0.5_dp*sqrt(2.0_dp*ev/(mi*Ti1))*dTi1ds + qi = Z1_global * qe + mi = am1_global * mu + vth = sqrt(2.0_dp * Ti1 * ev / mi) + dvthds = 0.5_dp * sqrt(2.0_dp * ev / (mi * Ti1)) * dTi1ds ! Call collision routine v0 = vth amb = 2.0_dp Zb = 1.0_dp - ebeam = amb*pmass*v0**2/(2.0_dp*ev) + ebeam = amb * pmass * v0**2 / (2.0_dp * ev) call loacol_nbi(amb, am1_global, am2_global, Zb, Z1_global, Z2_global, & ni1, ni2, Ti1, Ti2, Te, ebeam, v0, & dchichi, slowrate, dchichi_norm, slowrate_norm) @@ -210,8 +210,9 @@ subroutine init_thermodynamic_forces(psi_pr, q) real(dp), intent(in) :: psi_pr ! toroidal flux at plasma boundary == dpsi_tor/ds real(dp), intent(in) :: q ! safety factor - A1 = dni1ds/ni1 - qi/(Ti1*ev)*sign_theta*psi_pr/(q*c)*Om_tE - 3.0_dp/2.0_dp*dTi1ds/Ti1 - A2 = dTi1ds/Ti1 + A1 = dni1ds / ni1 - qi / (Ti1 * ev) * sign_theta * psi_pr / (q * c) * Om_tE - 3.0_dp / & + 2.0_dp * dTi1ds / Ti1 + A2 = dTi1ds / Ti1 end subroutine init_thermodynamic_forces end module neort_profiles diff --git a/src/resonance.f90 b/src/resonance.f90 index 4603bf39..143c336e 100644 --- a/src/resonance.f90 +++ b/src/resonance.f90 @@ -19,16 +19,16 @@ subroutine driftorbit_coarse(v, eta_min, eta_max, roots, nroots) ninterv = size(roots, 1) - deta = (eta_max - eta_min)*1.0_dp/ninterv + deta = (eta_max - eta_min) * 1.0_dp / ninterv nroots = 0 do k = 0, ninterv - eta = eta_min + k*deta + eta = eta_min + k * deta call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) call Om_ph(v, eta, Omph, dOmphdv, dOmphdeta) - res = mth*Omth + mph*Omph - dresdv = mth*dOmthdv + mph*dOmphdv - dresdeta = mth*dOmthdeta + mph*dOmphdeta + res = mth * Omth + mph * Omph + dresdv = mth * dOmthdv + mph * dOmphdv + dresdeta = mth * dOmthdeta + mph * dOmphdeta if (k > 0) then if (sign(1.0_dp, res) /= sign(1.0_dp, resold)) then nroots = nroots + 1 @@ -106,27 +106,27 @@ function driftorbit_root(v, tol, eta_min, eta_max) res_old = res call Om_ph(v, eta, Omph, dOmphdv, dOmphdeta) call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) - res = mph*Omph + mth*Omth + res = mph * Omph + mth * Omth driftorbit_root(1) = eta if (abs(res) < tol) then state = 1 - driftorbit_root(2) = mph*dOmphdeta + mth*dOmthdeta + driftorbit_root(2) = mph * dOmphdeta + mth * dOmthdeta exit elseif ((slope_pos .and. res > 0) .or. & ((.not. slope_pos) .and. res < 0)) then etamax2 = eta eta_old = eta - eta = (eta + etamin2)/2.0_dp + eta = (eta + etamin2) / 2.0_dp else etamin2 = eta eta_old = eta - eta = (eta + etamax2)/2.0_dp + eta = (eta + etamax2) / 2.0_dp end if end do if (state < 0) then - driftorbit_root(2) = mph*dOmphdeta + mth*dOmthdeta + driftorbit_root(2) = mph * dOmphdeta + mth * dOmthdeta write (msg, "(a,i0,a,g0,a,i0,a,g0,a,g0,a,g0,a,g0,a,g0,a,g0,a,g0,a,g0,a,g0,a,g0)") & "driftorbit_root: did not converge within ", maxit, " iterations" // LF // & TAB // "v/vth = ", v / vth, ", mth = ", mth, ", sign_vpar = ", sign_vpar, LF // & diff --git a/src/transport.f90 b/src/transport.f90 index 92599bb3..fe74171b 100644 --- a/src/transport.f90 +++ b/src/transport.f90 @@ -49,31 +49,32 @@ pure function D11int(ux, taub, Hmn2) real(dp) :: D11int real(dp), intent(in) :: ux, taub, Hmn2 - D11int = pi**(3.0_dp/2.0_dp)*mph**2*c**2*q*vth & - /(qi**2*dVds*abs(psi_pr))*ux**3*exp(-ux**2)*taub*Hmn2 + D11int = pi**(3.0_dp / 2.0_dp) * mph**2 * c**2 * q * vth & + / (qi**2 * dVds * abs(psi_pr)) * ux**3 * exp(-ux**2) * taub * Hmn2 end function D11int pure function D12int(ux, taub, Hmn2) real(dp) :: D12int real(dp), intent(in) :: ux, taub, Hmn2 - D12int = D11int(ux, taub, Hmn2)*ux**2 + D12int = D11int(ux, taub, Hmn2) * ux**2 end function D12int pure function Tphi_int(ux, taub, Hmn2) real(dp) :: Tphi_int real(dp), intent(in) :: ux, taub, Hmn2 - Tphi_int = sign(1.0_dp, psi_pr*q*sign_theta)*pi**(3.0_dp/2.0_dp)*mph**2*ni1*c*vth/qi & - *ux**3*exp(-ux**2)*taub*Hmn2*(A1 + A2*ux**2) + Tphi_int = sign(1.0_dp, psi_pr * q * sign_theta) * pi**(3.0_dp / 2.0_dp) * mph**2 * ni1 * & + c * vth / qi & + * ux**3 * exp(-ux**2) * taub * Hmn2 * (A1 + A2 * ux**2) end function Tphi_int subroutine compute_transport_integral(vmin, vmax, vsteps, D, T) ! compute transport integral via midpoint rule real(dp), intent(in) :: vmin, vmax integer, intent(in) :: vsteps - real(dp), intent(out) :: D(2), T ! Transport coefficients D and torque density T - real(dp) :: D_plateau, dsdreff ! Plateau diffusion coefficient and ds/dreff=<|grad s|> + real(dp), intent(out) :: D(2), T ! Transport coefficients D and torque density T + real(dp) :: D_plateau, dsdreff ! Plateau diffusion coefficient and ds/dreff=<|grad s|> real(dp) :: ux, du, dD11, dD12, dT, v, eta real(dp) :: eta_res(2) real(dp) :: taub, bounceavg(nvar) @@ -86,20 +87,20 @@ subroutine compute_transport_integral(vmin, vmax, vsteps, D, T) D = 0.0_dp T = 0.0_dp - du = (vmax - vmin)/(vsteps*vth) - ux = vmin/vth + du/2.0_dp + du = (vmax - vmin) / (vsteps * vth) + ux = vmin / vth + du / 2.0_dp do ku = 1, vsteps - v = ux*vth + v = ux * vth call driftorbit_coarse(v, etamin, etamax, roots, nroots) if (nroots == 0) continue do kr = 1, nroots - eta_res = driftorbit_root(v, 1.0e-8_dp*abs(Om_tE), roots(kr, 1), roots(kr, 2)) + eta_res = driftorbit_root(v, 1.0e-8_dp * abs(Om_tE), roots(kr, 1), roots(kr, 2)) eta = eta_res(1) call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) - taub = 2.0_dp*pi/abs(Omth) + taub = 2.0_dp * pi / abs(Omth) call bounce_fast(v, eta, taub, bounceavg, timestep_transport, istate_dv) if (istate_dv == -1) then call error(fmt_dbg('VODE MXSTEP: mth=', dble(mth), ' ux=', ux, ' eta=', eta, ' taub=', taub)) @@ -110,26 +111,26 @@ subroutine compute_transport_integral(vmin, vmax, vsteps, D, T) call trace(fmt_dbg('near etatp: mth=', dble(mth), ' ux=', ux, ' eta=', eta, ' taub=', taub)) end if end if - Hmn2 = (bounceavg(3)**2 + bounceavg(4)**2)*(mi*(ux*vth)**2/2.0_dp)**2 + Hmn2 = (bounceavg(3)**2 + bounceavg(4)**2) * (mi * (ux * vth)**2 / 2.0_dp)**2 attenuation_factor = nonlinear_attenuation(ux, eta, bounceavg, Omth, & dOmthdv, dOmthdeta, Hmn2) - dD11 = du*D11int(ux, taub, Hmn2)/abs(eta_res(2)) - dD12 = du*D12int(ux, taub, Hmn2)/abs(eta_res(2)) - D(1) = D(1) + dD11*attenuation_factor - D(2) = D(2) + dD12*attenuation_factor + dD11 = du * D11int(ux, taub, Hmn2) / abs(eta_res(2)) + dD12 = du * D12int(ux, taub, Hmn2) / abs(eta_res(2)) + D(1) = D(1) + dD11 * attenuation_factor + D(2) = D(2) + dD12 * attenuation_factor if (comptorque) then - dT = du*Tphi_int(ux, taub, Hmn2)/abs(eta_res(2)) - T = T + dT*attenuation_factor + dT = du * Tphi_int(ux, taub, Hmn2) / abs(eta_res(2)) + T = T + dT * attenuation_factor end if end do ux = ux + du end do - D_plateau = pi*vth**3/(16.0_dp*R0*iota*(qi*B0/(mi*c))**2) - dsdreff = 2.0_dp/a*sqrt(s) ! TODO: Use exact value instead of this approximation - D = dsdreff**(-2)*D/D_plateau + D_plateau = pi * vth**3 / (16.0_dp * R0 * iota * (qi * B0 / (mi * c))**2) + dsdreff = 2.0_dp / a * sqrt(s) ! TODO: Use exact value instead of this approximation + D = dsdreff**(-2) * D / D_plateau call debug("compute_transport_integral complete") @@ -152,7 +153,7 @@ subroutine timestep_transport(v, eta, neq, t, y, ydot) ! for y(1:3) real(dp) :: bmod, sqrtg, x(3), hder(3), hcovar(3), hctrvr(3), hcurl(3), Om_tB_v real(dp) :: t0 - complex(dp) :: epsn, Hn ! relative amplitude of perturbation field epsn=Bn/B0 + complex(dp) :: epsn, Hn ! relative amplitude of perturbation field epsn=Bn/B0 ! and Hamiltonian Hn = (H - H0)_n x(1) = s @@ -164,24 +165,26 @@ subroutine timestep_transport(v, eta, neq, t, y, ydot) ! evaluate orbit averages of Hamiltonian perturbation if (pertfile) then call do_magfie_pert_amp(x, epsn) - epsn = epsmn*epsn/bmod + epsn = epsmn * epsn / bmod else - epsn = epsmn*exp(imun*m0*y(1)) + epsn = epsmn * exp(imun * m0 * y(1)) end if if (eta > etatp) then !t0 = 0.25*2*pi/Omth ! Different starting position in orbit t0 = 0.0_dp - Hn = (2.0_dp - eta*bmod)*epsn*exp(imun*(q*mph*(y(1)) - mth*(t - t0)*Omth)) + Hn = (2.0_dp - eta * bmod) * epsn * exp(imun * (q * mph * (y(1)) - mth * (t - & + t0) * Omth)) else - Hn = (2.0_dp - eta*bmod)*epsn*exp(imun*(q*mph*(y(1)) - (mth + q*mph)*t*Omth)) + Hn = (2.0_dp - eta * bmod) * epsn * exp(imun * (q * mph * (y(1)) - (mth + q * mph) & + * t * Omth)) end if ydot(3) = real(Hn) ydot(4) = aimag(Hn) ! evaluate orbit averages for nonlinear attenuation if (nonlin) then - ydot(5) = 1.0_dp/bmod + ydot(5) = 1.0_dp / bmod ydot(6) = bmod else ydot(5:6) = 0.0_dp diff --git a/src/util.f90 b/src/util.f90 index 8403c0ce..0767c587 100644 --- a/src/util.f90 +++ b/src/util.f90 @@ -10,17 +10,17 @@ module util complex(dp), parameter :: imun = (0.0_dp, 1.0_dp) - real(dp), parameter, public :: pi = 4*atan(1.0_dp) + real(dp), parameter, public :: pi = 4 * atan(1.0_dp) real(dp), parameter, public :: & - qe = 4.803204e-10_dp, & ! elementary charge + qe = 4.803204e-10_dp, & ! elementary charge !me = 9.109382e-28_dp, & ! electron mass, - mu = 1.660538e-24_dp, & ! 1u - c = 2.997925e+10_dp, & ! speed of light + mu = 1.660538e-24_dp, & ! 1u + c = 2.997925e+10_dp, & ! speed of light !kb = 1.381649e-16_dp, & ! Boltzmann constant - eV = 1.602176e-12_dp ! 1 electron volt + eV = 1.602176e-12_dp ! 1 electron volt - real(dp), public :: qi = 1.0_dp*qe, mi = 2.014_dp*mu + real(dp), public :: qi = 1.0_dp * qe, mi = 2.014_dp * mu ! Flux-surface dependent (set from plasma data at each s) !$omp threadprivate (qi, mi) diff --git a/test/test_bounce.f90 b/test/test_bounce.f90 index d69dd97e..7dd9ec4a 100644 --- a/test/test_bounce.f90 +++ b/test/test_bounce.f90 @@ -23,15 +23,15 @@ subroutine setup call do_magfie_init("in_file") call init - Om_tE = vth*M_t/R0 - dOm_tEds = vth*dM_tds/R0 + Om_tE = vth * M_t / R0 + dOm_tEds = vth * dM_tds / R0 - etamin = (1+epst)*etatp - etamax = (1-epst)*etadt + etamin = (1 + epst) * etatp + etamax = (1 - epst) * etadt sign_vpar = 1 v = vth - eta = 0.5_dp*(etamin + etamax) + eta = 0.5_dp * (etamin + etamax) end subroutine setup diff --git a/test/test_misc.f90 b/test/test_misc.f90 index cde8482b..f66f96ff 100644 --- a/test/test_misc.f90 +++ b/test/test_misc.f90 @@ -35,10 +35,10 @@ subroutine test_torfreq v = vth etamin = etatp - etamax = etatp + (etadt - etatp)*(1.0_dp - epsst_spl) + etamax = etatp + (etadt - etatp) * (1.0_dp - epsst_spl) b = log(epst_spl) - aa = 1.0_dp/(n - 1.0_dp)*(log(etamax/etamin - 1.0_dp) - b) + aa = 1.0_dp / (n - 1.0_dp) * (log(etamax / etamin - 1.0_dp) - b) eta = etamax sign_vpar = 1 @@ -98,8 +98,8 @@ subroutine test_resline print *, "test_resline" - vmin = 1.0e-6_dp*vth - vmax = 10.0_dp*vth + vmin = 1.0e-6_dp * vth + vmax = 10.0_dp * vth etaresp = etatp etarest = etatp @@ -109,7 +109,7 @@ subroutine test_resline open (unit=10, file=trim(adjustl(runname))//"_resline_pct.out", recl=1024) open (unit=11, file=trim(adjustl(runname))//"_resline_t.out", recl=1024) do k = 0, n - 1 - v = vmin + k/(n - 1.0_dp)*(vmax - vmin) + v = vmin + k / (n - 1.0_dp) * (vmax - vmin) if (.not. nopassing) then ! resonance (passing) diff --git a/test/test_omega_prime.f90 b/test/test_omega_prime.f90 index 55fd29bc..9329e1fb 100644 --- a/test/test_omega_prime.f90 +++ b/test/test_omega_prime.f90 @@ -45,7 +45,7 @@ subroutine main call read_and_set_config(runname) call setup s0 = s - v0 = 0.9_dp*vth + v0 = 0.9_dp * vth eta_res = first_resonance(v0) eta0 = eta_res(1) !eta0 = 0.5_dp*(etamin + etamax) @@ -62,19 +62,19 @@ subroutine sample_values(s0, v0, eta0, freq_data) integer :: i real(dp) :: xi - freq_data(1) % s = s0 - s = freq_data(1) % s - freq_data(1) % v = v0 - freq_data(1) % eta = eta0 + freq_data(1)%s = s0 + s = freq_data(1)%s + freq_data(1)%v = v0 + freq_data(1)%eta = eta0 call test_omega_prime(freq_data(1)) do i = 2, size(freq_data) call random_number(xi) - freq_data(i) % s = s0*(1.0_dp + DELTA*(xi - 0.5_dp)) - s = freq_data(i) % s + freq_data(i)%s = s0 * (1.0_dp + DELTA * (xi - 0.5_dp)) + s = freq_data(i)%s call random_number(xi) - freq_data(i) % v = v0*(1.0_dp + DELTA*(xi - 0.5_dp)) + freq_data(i)%v = v0 * (1.0_dp + DELTA * (xi - 0.5_dp)) call random_number(xi) - freq_data(i) % eta = eta0 + (etamax - etamin)*DELTA*(xi - 0.5_dp) + freq_data(i)%eta = eta0 + (etamax - etamin) * DELTA * (xi - 0.5_dp) call test_omega_prime(freq_data(i)) end do end subroutine sample_values @@ -124,7 +124,7 @@ function first_resonance(v) result(eta_res) real(dp) :: roots(nlev, 3) integer :: nroots call driftorbit_coarse(v, etamin, etamax, roots, nroots) - eta_res = driftorbit_root(v, 1.0e-8_dp*abs(Om_tE), roots(1, 1), roots(1, 2)) + eta_res = driftorbit_root(v, 1.0e-8_dp * abs(Om_tE), roots(1, 1), roots(1, 2)) end function first_resonance subroutine test_omega_prime(f) @@ -135,19 +135,19 @@ subroutine test_omega_prime(f) call compute_frequencies(f) - print *, "dOmthds_test: ", dOmthds_test(f % v, f % eta) + print *, "dOmthds_test: ", dOmthds_test(f%v, f%eta) - call bounce_fast(f % v, f % eta, 2.0_dp*pi/abs(f % Omth), bounceavg, timestep_transport) + call bounce_fast(f%v, f%eta, 2.0_dp * pi / abs(f%Omth), bounceavg, timestep_transport) - f % Ompr_old = omega_prime_old(f % v/vth, f % eta, f % Omth, f % dOmdv, f % dOmdeta, f % dOmdpph) - f % Ompr_new = omega_prime(f % v/vth, f % eta, bounceavg, f % Omth, f % dOmdv, f % dOmdeta, f % dOmdpph) + f%Ompr_old = omega_prime_old(f%v / vth, f%eta, f%Omth, f%dOmdv, f%dOmdeta, f%dOmdpph) + f%Ompr_new = omega_prime(f%v / vth, f%eta, bounceavg, f%Omth, f%dOmdv, f%dOmdeta, f%dOmdpph) - call compute_invariants(f % v, f % eta, f % J) + call compute_invariants(f%v, f%eta, f%J) ! Transformed actions for Galileo transformation near resonance - f % Jbar(1) = f % J(1) - f % Jbar(2) = f % J(2) - mth/mph*f % J(3) - f % Jbar(3) = f % J(3)/mph + f%Jbar(1) = f%J(1) + f%Jbar(2) = f%J(2) - mth / mph * f%J(3) + f%Jbar(3) = f%J(3) / mph end subroutine test_omega_prime subroutine compute_frequencies(f) @@ -155,18 +155,18 @@ subroutine compute_frequencies(f) real(dp) :: taub, bounceavg(nvar) - call Om_th(f % v, f % eta, f % Omth, f % dOmthdv, f % dOmthdeta) - taub = 2.0_dp*pi/abs(f % Omth) - call bounce_fast(f % v, f % eta, taub, bounceavg, timestep_transport) - call Om_ph(f % v, f % eta, f % Omph, f % dOmphdv, f % dOmphdeta) - call d_Om_ds(f % v, f % eta, taub, f % dOmthds, f % dOmphds) - f % Om = mth*f % Omth + mph*f % Omph - f % dOmdv = mth*f % dOmthdv + mph*f % dOmphdv - f % dOmdeta = mth*f % dOmthdeta + mph*f % dOmphdeta - f % dOmds = mth*f % dOmthds + mph*f % dOmphds - f % dOmdpph = -(qi/c*iota*sign_theta*psi_pr)**(-1)*f % dOmds - f % Om_tE = Om_tE - f % dOm_tEds = dOm_tEds + call Om_th(f%v, f%eta, f%Omth, f%dOmthdv, f%dOmthdeta) + taub = 2.0_dp * pi / abs(f%Omth) + call bounce_fast(f%v, f%eta, taub, bounceavg, timestep_transport) + call Om_ph(f%v, f%eta, f%Omph, f%dOmphdv, f%dOmphdeta) + call d_Om_ds(f%v, f%eta, taub, f%dOmthds, f%dOmphds) + f%Om = mth * f%Omth + mph * f%Omph + f%dOmdv = mth * f%dOmthdv + mph * f%dOmphdv + f%dOmdeta = mth * f%dOmthdeta + mph * f%dOmphdeta + f%dOmds = mth * f%dOmthds + mph * f%dOmphds + f%dOmdpph = -(qi / c * iota * sign_theta * psi_pr)**(-1) * f%dOmds + f%Om_tE = Om_tE + f%dOm_tEds = dOm_tEds end subroutine compute_frequencies subroutine compute_invariants(v, eta, J) @@ -183,13 +183,13 @@ subroutine compute_invariants(v, eta, J) ! is defined inside subroutine timestep (thin orbit integration) call evaluate_bfield_local(bmod, htheta) y0(1) = th0 - y0(2) = sign(1.0_dp, htheta)*sign_vpar*vpar(v, eta, bmod) + y0(2) = sign(1.0_dp, htheta) * sign_vpar * vpar(v, eta, bmod) y0(3) = 1.0e-15_dp - taub = 2.0_dp*pi/(vperp(v, eta, bmod)*iota/R0*sqrt(eps/2.0_dp)) - bounceint = bounce_integral(v, eta, neq, y0, taub/5.0_dp, timestep_invariants) + taub = 2.0_dp * pi / (vperp(v, eta, bmod) * iota / R0 * sqrt(eps / 2.0_dp)) + bounceint = bounce_integral(v, eta, neq, y0, taub / 5.0_dp, timestep_invariants) taub = bounceint(1) - bounceavg = bounceint(2:)/taub + bounceavg = bounceint(2:) / taub print *, "taub: ", taub print *, "bounceavg: ", bounceavg @@ -219,18 +219,18 @@ subroutine timestep_invariants(v, eta, neq, t, y, ydot) call do_magfie(x, bmod, sqrtg, hder, hcovar, hctrvr, hcurl) call poloidal_velocity(v, eta, bmod, hctrvr(3), hder(3), y(2), ydot(:2)) - ydot(3) = mi*ydot(2)**2 ! Jpar + ydot(3) = mi * ydot(2)**2 ! Jpar end subroutine timestep_invariants pure function Jperp(v, eta) real(dp) :: Jperp real(dp), intent(in) :: v, eta - Jperp = 0.5_dp*mi*v**2*mi*c/qi*eta + Jperp = 0.5_dp * mi * v**2 * mi * c / qi * eta end function Jperp function pphi() real(dp) :: pphi - pphi = -qi/c*psi_pol() + pphi = -qi / c * psi_pol() end function pphi function psi_pol() @@ -279,7 +279,7 @@ function trapz_step(x1, x2, y1, y2) real(dp), intent(in) :: x1, x2, y1, y2 real(dp) :: trapz_step - trapz_step = 0.5_dp*(x2 - x1)*(y1 + y2) + trapz_step = 0.5_dp * (x2 - x1) * (y1 + y2) end function trapz_step function dOmthds_test(v, eta) @@ -293,12 +293,12 @@ function dOmthds_test(v, eta) s0 = s s = s0 + ds - Omth1 = 2.0_dp*pi/bounce_time(v, eta) + Omth1 = 2.0_dp * pi / bounce_time(v, eta) s = s0 - ds - Omth2 = 2.0_dp*pi/bounce_time(v, eta) + Omth2 = 2.0_dp * pi / bounce_time(v, eta) - dOmthds_test = (Omth1 - Omth2)/(2.0_dp*ds) + dOmthds_test = (Omth1 - Omth2) / (2.0_dp * ds) s = s0 end function dOmthds_test @@ -306,7 +306,8 @@ end function dOmthds_test function omega_prime_old(ux, eta, Omth, dOmdv, dOmdeta, dOmdpph) real(dp), intent(in) :: ux, eta, Omth, dOmdv, dOmdeta, dOmdpph real(dp) :: omega_prime_old - omega_prime_old = mth*(eta*dOmdeta - ux*vth/2*dOmdv)/(mi*(ux*vth)**2/(2.0_dp*Omth)) + dOmdpph + omega_prime_old = mth * (eta * dOmdeta - ux * vth / 2 * dOmdv) / (mi * (ux * vth)**2 / & + (2.0_dp * Omth)) + dOmdpph end function omega_prime_old end program test_omage_prime_prog diff --git a/test/test_torque.f90 b/test/test_torque.f90 index e28bd23e..f597fc37 100644 --- a/test/test_torque.f90 +++ b/test/test_torque.f90 @@ -17,7 +17,7 @@ program test_torque_prog integer, parameter :: MTH = -3 real(dp), parameter :: AVG_NABLA_S = 0.02162413513580247_dp real(dp), parameter :: DV_OVER_DS = 12621249.179676009_dp - real(dp), parameter :: FLUX_SURFACE_AREA = DV_OVER_DS*AVG_NABLA_S + real(dp), parameter :: FLUX_SURFACE_AREA = DV_OVER_DS * AVG_NABLA_S real(dp), parameter :: TOL = 1.0e-7_dp character(1024) :: TEST_RUN = "driftorbit64.new" type(magfie_data_t) :: magfie_data @@ -66,14 +66,14 @@ subroutine test_torque write (*, "(4ES12.2,2F12.2)") TTco, TTctr, TTt, TTco + TTctr + TTt print *, "" print *, "Ratios to direct calculation:" - print *, "TTctr/Tctr = ", TTctr/Tctr - print *, "TTt/Tt = ", TTt/Tt + print *, "TTctr/Tctr = ", TTctr / Tctr + print *, "TTt/Tt = ", TTt / Tt - if (abs(TTctr/Tctr - 1.0_dp) > TOL) then + if (abs(TTctr / Tctr - 1.0_dp) > TOL) then error stop end if - if (abs(TTt/Tt - 1.0_dp) > TOL) then + if (abs(TTt / Tt - 1.0_dp) > TOL) then error stop end if @@ -92,11 +92,11 @@ function torque_from_transport(D) result(Tphi) x(2) = 0.0_dp x(3) = 0.0_dp call do_magfie(x, bmod, sqrtg, hder, hcovar, hctrvr, hcurl) - sqrtgBth = sqrtg*hctrvr(3)*bmod*AVG_NABLA_S + sqrtgBth = sqrtg * hctrvr(3) * bmod * AVG_NABLA_S print *, "sqrtgBth = ", sqrtgBth - D_plateau = pi*vth**3/(16.0_dp*R0*iota*(qi*B0/(mi*c))**2) - Tphi = (sqrtgBth/c)*qe*(-ni1*D_plateau*(D(1)*A1 + D(2)*A2))*AVG_NABLA_S + D_plateau = pi * vth**3 / (16.0_dp * R0 * iota * (qi * B0 / (mi * c))**2) + Tphi = (sqrtgBth / c) * qe * (-ni1 * D_plateau * (D(1) * A1 + D(2) * A2)) * AVG_NABLA_S Tphi = Tphi end function torque_from_transport @@ -104,13 +104,12 @@ end function torque_from_transport ! This then finally cancels out in torque calculation subroutine correct_transport_coeff(D) real(dp), intent(inout) :: D(2) - D = D*(2.0_dp/a*sqrt(s))**2/AVG_NABLA_S**2 + D = D * (2.0_dp / a * sqrt(s))**2 / AVG_NABLA_S**2 end subroutine correct_transport_coeff - ! Correct torque calculation subroutine correct_torque(T) real(dp), intent(inout) :: T - T = T/DV_OVER_DS + T = T / DV_OVER_DS end subroutine correct_torque end program test_torque_prog From 2e121a62ad1e21f9546e2c6f2a587446085a9755 Mon Sep 17 00:00:00 2001 From: Andreas Zach <80193776+zandivx@users.noreply.github.com> Date: Sat, 10 Jan 2026 01:35:31 +0100 Subject: [PATCH 3/3] Update test/test_reslines.f90 Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- test/test_reslines.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_reslines.f90 b/test/test_reslines.f90 index 18151cdc..2b69a505 100644 --- a/test/test_reslines.f90 +++ b/test/test_reslines.f90 @@ -10,7 +10,7 @@ program test_reslines use neort_freq, only: Om_th, Om_ph implicit none - real(dp), parameter :: vnorm = sqrt(0.15d0) + real(dp), parameter :: vnorm = sqrt(0.15_dp) integer :: fid