Skip to content

Make constants explicitly double precision, otherwise they are unintentionally treated as single precision by Fortran 90 #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Make constants explicitly double precision, otherwise they are unintentionally treated as single precision by Fortran 90 #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
22e827c
Made double precision explicit in constants
kryl0v-subspace Jul 15, 2021
fb2387e
Explicit _dp tweaks
kryl0v-subspace Jul 19, 2021
d6c4dd2
Made pi explicitly double precision
kryl0v-subspace Jul 19, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 12 additions & 12 deletions c2ray_parameters.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ module c2ray_parameters

!> Which fraction of the cells can be left unconverged in order
!! to improve performance (used in rad_evolve3d)
real(kind=dp),parameter :: convergence_fraction=1.5e-5
real(kind=dp),parameter :: convergence_fraction=1.5e-5_dp

! Set to true to let C2-Ray not change the temperature
! set interactively in mat_ini
Expand All @@ -31,13 +31,13 @@ module c2ray_parameters
real(kind=dp),parameter :: epsilon=1e-40_dp

!> Convergence criterion for per source calculation (evolve0d)
real(kind=dp),parameter :: convergence1=1.0e-3
real(kind=dp),parameter :: convergence1=1.0e-3_dp

!> Convergence criterion for global calculation (evolve0d)
real(kind=dp),parameter :: convergence2=1.0e-2
real(kind=dp),parameter :: convergence2=1.0e-2_dp

!> Parameters for nominal SED
real(kind=dp),parameter :: teff_nominal=0.0
real(kind=dp),parameter :: teff_nominal=0.0_dp
real(kind=dp),parameter :: s_star_nominal=1e48_dp
!real(kind=dp),parameter :: s_star_nominal=1e50_dp

Expand All @@ -48,27 +48,27 @@ module c2ray_parameters
!! 4: 1 Mpc P3M
integer,parameter :: type_of_clumping=3
!> Clumping factor if constant
real,parameter :: clumping_factor=1.0
real,parameter :: clumping_factor=1.0_dp

! Cosmological cooling
! Set in cosmology module!!
!logical,parameter :: cosmological=.true.

!> Thermal: minimum temperature
real(kind=dp),parameter :: minitemp=1.0 ! minimum temperature
real(kind=dp),parameter :: minitemp=1.0_dp ! minimum temperature
!> Thermal: fraction of the cooling time step below which no iteration is done
real(kind=dp),parameter :: relative_denergy=0.1
real(kind=dp),parameter :: relative_denergy=0.1_dp

!
!> Source properties: Photon per atom for high mass sources
real,parameter :: phot_per_atom1=250.0
real,parameter :: phot_per_atom1=250.0_dp
!> Source properties: Photon per atom for low mass sources
real,parameter :: phot_per_atom2=250.0
real,parameter :: phot_per_atom2=250.0_dp
!> Source properties: Life time of sources (if set at compile time)
real,parameter :: lifetime=20e6*YEAR
real,parameter :: lifetime=20e6_dp*YEAR
!> Source properties: Lower limit neutral fraction for suppression criterion
real,parameter :: StillNeutral=0.9 ! lower limit of neutral criterium
real,parameter :: StillNeutral=0.9_dp ! lower limit of neutral criterium
!> Source properties: Upper limit for low mass sources (not used)
real,parameter :: LowMassLimit=1e9 ! in solar masses
real,parameter :: LowMassLimit=1e9_dp ! in solar masses

end module c2ray_parameters
16 changes: 8 additions & 8 deletions cgsastroconstants.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,16 +20,16 @@ module astroconstants
! A collection of astronomical units and conversion factors
! Units: cgs

real(kind=dp),parameter :: R_SOLAR=6.9599e10 !< Solar radius
real(kind=dp),parameter :: L_SOLAR=3.826e33 !< Solar luminosity
real(kind=dp),parameter :: R_SOLAR=6.9599e10_dp !< Solar radius
real(kind=dp),parameter :: L_SOLAR=3.826e33_dp !< Solar luminosity
real(kind=dp),parameter :: M_SOLAR=1.98892d33 !< Solar mass

real(kind=dp),parameter :: YEAR=3.15576E+07 !< Julian year
real(kind=dp),parameter :: YEAR=3.15576E+07_dp !< Julian year

real(kind=dp),parameter :: pc=3.086e18 !< parsec
real(kind=dp),parameter :: kpc=1e3*pc !< kiloparsec
real(kind=dp),parameter :: Mpc=1e6*pc !< megaparsec
real(kind=dp),parameter :: lightyear=9.463e17 !< lightyear
real(kind=dp),parameter :: AU=1.49597870E+13 !< Astronomical Unit
real(kind=dp),parameter :: pc=3.086e18_dp !< parsec
real(kind=dp),parameter :: kpc=1e3_dp*pc !< kiloparsec
real(kind=dp),parameter :: Mpc=1e6_dp*pc !< megaparsec
real(kind=dp),parameter :: lightyear=9.463e17_dp !< lightyear
real(kind=dp),parameter :: AU=1.49597870E+13_dp !< Astronomical Unit

end module astroconstants
32 changes: 16 additions & 16 deletions cgsconstants.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,31 +34,31 @@ module cgsconstants
real(kind=dp), parameter :: G_grav=6.6732e-8_dp

!> ev2k - conversion factor between evs and kelvins
real(kind=dp),parameter :: ev2k=1.0/8.617e-05
real(kind=dp),parameter :: ev2k=1.0_dp/8.617e-05_dp
!> ev2erg - conversion factor between evs and ergs
real(kind=dp),parameter :: ev2erg=1.602e-12
real(kind=dp),parameter :: ev2erg=1.602e-12_dp
!> ev2j - conversion factor between ergs and Joules
real(kind=dp),parameter :: erg2j=1e-7
real(kind=dp),parameter :: erg2j=1e-7_dp

! The following are scaled to frequency scaling

!> Frequency scaling factor,
!! this scaling parameter is independent of any main program scaling
!! (see scaling.f90), and may only be used in the radiation physics
!! subroutines (currently switched off)
real(kind=dp),parameter :: sclfre=1.0e15
real(kind=dp),parameter :: sclfre=1.0e15_dp
!> conversion between evs and frequency
real(kind=dp),parameter :: ev2fr=0.241838e15!/sclfre
real(kind=dp),parameter :: ev2fr=0.241838e15_dp!/sclfre

! h/k, Planck/Boltzm
! Check this number, seems scaled
!real(kind=dp),parameter :: hoverk=47979.72484
!> Planck constant scaled
real(kind=dp),parameter :: hscl=hplanck!*sclfre
!> tpic2 - 2*pi/c^2 times scaling factors needed for the integral cores
real(kind=dp),parameter :: tpic2=2.0*pi/(c*c)
real(kind=dp),parameter :: tpic2=2.0_dp*pi/(c*c)
!> two_pi_c2 - 2*pi/c^2 times scaling factors needed for the integral cores
real(kind=dp),parameter :: two_pi_c2=2.0*pi/(c*c)!*sclfre**3
real(kind=dp),parameter :: two_pi_c2=2.0_dp*pi/(c*c)!*sclfre**3

!> Hydrogen recombination parameter (power law index)
real(kind=dp),parameter :: albpow=-0.7_dp !in the old code -0.79
Expand All @@ -73,32 +73,32 @@ module cgsconstants
!here it was 5.91e-12 I replace with book value of 2*1.1e-12

!> Hydrogen ionization energy (in eV)
real(kind=dp), parameter :: eth0=13.598
real(kind=dp), parameter :: eth0=13.598_dp
!> Hydrogen ionization energy (in erg)
real(kind=dp),parameter :: hionen=eth0*ev2erg
!> Hydrogen ionization energy expressed in K
real(kind=dp),parameter :: temph0=eth0*ev2k
!> Hydrogen collisional ionization parameter 1
real(kind=dp),parameter :: xih0=1.0
real(kind=dp),parameter :: xih0=1.0_dp
!> Hydrogen collisional ionization parameter 2
real(kind=dp),parameter :: fh0=0.83
real(kind=dp),parameter :: fh0=0.83_dp
!> Hydrogen collisional ionization parameter
real(kind=dp),parameter :: colh0=1.3e-8*fh0*xih0/(eth0*eth0)
real(kind=dp),parameter :: colh0=1.3e-8_dp*fh0*xih0/(eth0*eth0)


!> Helium ionization energy (in eV)
real(kind=dp), dimension(0:1), parameter :: ethe=(/24.587,54.416/)
real(kind=dp), dimension(0:1), parameter :: ethe=(/24.587_dp,54.416_dp/)
!> Helium ionization energy (in erg)
real(kind=dp), dimension(0:1), parameter :: heionen=(/ethe(0)*ev2erg,ethe(1)*ev2erg/)
!> Helium ionization energy expressed in K
real(kind=dp), dimension(0:1), parameter :: temphe=(/ethe(0)*ev2k,ethe(1)*ev2k/)
!> Helium collisional ionization parameter 1
real(kind=dp),dimension(0:1),parameter :: xihe=(/2.0,1.0/)
real(kind=dp),dimension(0:1),parameter :: xihe=(/2.0_dp,1.0_dp/)
!> Helium collisional ionization parameter 2
real(kind=dp),dimension(0:1),parameter :: fhe=(/0.63,1.30/)
real(kind=dp),dimension(0:1),parameter :: fhe=(/0.63_dp,1.30_dp/)
!> Helium collisional ionization parameter
real(kind=dp),dimension(0:1),parameter :: colhe=(/1.3e-8*fhe(0)*xihe(0)/(ethe(0)*ethe(0)), &
1.3e-8*fhe(1)*xihe(1)/(ethe(1)*ethe(1))/)
real(kind=dp),dimension(0:1),parameter :: colhe=(/1.3e-8_dp*fhe(0)*xihe(0)/(ethe(0)*ethe(0)), &
1.3e-8_dp*fhe(1)*xihe(1)/(ethe(1)*ethe(1))/)



Expand Down
24 changes: 12 additions & 12 deletions cgsphotoconstants.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,41 +22,41 @@ module cgsphotoconstants
! !> Helium ionization potentials (eV)
! real(kind=dp), dimension(0:1),parameter :: ethe=(/24.587,54.416/)
!> Hydrogen cross section
real(kind=dp), parameter :: sigh=6.30e-18
real(kind=dp), parameter :: sigh=6.30e-18_dp
!> Helium cross section
real(kind=dp), parameter :: sighe0=7.83e-18
real(kind=dp), parameter :: sighe0=7.83e-18_dp
!> He+ cross section
real(kind=dp), parameter :: sighe1=1.58e-18
real(kind=dp), parameter :: sighe1=1.58e-18_dp
!> H ionization energy in frequency
real(kind=dp), parameter :: frth0=ev2fr*eth0
!> He ionization energy in frequency
real(kind=dp), parameter :: frthe0=ev2fr*ethe(0)
!> He+ ionization energy in frequency
real(kind=dp), parameter :: frthe1=ev2fr*ethe(1)
!> Frequency dependence of H cross section parameter
real(kind=dp),parameter :: betah0=1.0
real(kind=dp),parameter :: betah0=1.0_dp
!> Frequency dependence of He cross section parameter
real(kind=dp), dimension(0:1),parameter :: betahe=(/1.0,1.0/)
real(kind=dp), dimension(0:1),parameter :: betahe=(/1.0_dp,1.0_dp/)

!> Frequency dependence of H cross section parameter
real(kind=dp), parameter :: sh0=2.8
real(kind=dp), parameter :: sh0=2.8_dp
!> Frequency dependence of He cross section parameter
real(kind=dp), parameter :: she0=1.7
real(kind=dp), parameter :: she0=1.7_dp
!> Frequency dependence of He+ cross section parameter
real(kind=dp), parameter :: she1=2.8
real(kind=dp), parameter :: she1=2.8_dp

!> maximum T_eff for black body
real(kind=dp), parameter :: thigh=200000.0
real(kind=dp), parameter :: thigh=200000.0_dp
!> minimum T_eff for black body
real(kind=dp), parameter :: tlow=2000.0
real(kind=dp), parameter :: tlow=2000.0_dp
!real(kind=dp), parameter :: frtop1=700.0*tlow/47979.72484*1e15

!> Upper limits for integrals: due to arithmetic precision,
!! exp(700) exceeds double precision limit
real(kind=dp), parameter :: frtop1=700.0*tlow*kb/hplanck
real(kind=dp), parameter :: frtop1=700.0_dp*tlow*kb/hplanck
!> Upper limits for integrals: due to the form of the planck
!! curve: take 10 times the frequency of maximum intensity
real(kind=dp), parameter :: frtop2=5.88e-05*thigh*1e15
real(kind=dp), parameter :: frtop2=5.88e-05_dp*thigh*1e15_dp

!> H optical depth fit parameter for frequency range 2
real(kind=dp) :: tf2h
Expand Down
2 changes: 1 addition & 1 deletion mathconstants.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,6 @@ module mathconstants
private

!> the number pi
real(kind=dp),public,parameter :: pi=3.141592654
real(kind=dp),public,parameter :: pi=3.141592654_dp

end module mathconstants
4 changes: 2 additions & 2 deletions radiation.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -529,7 +529,7 @@ subroutine photoion (phi,hcolum_in,hcolum_out,vol)
! find the table positions for the optical depth (ingoing)
tau1=log10(max(1.0e-20_dp,tauh_in))
! odpos1=min(1.0d0*NumTau,max(0.0d0,1.0d0+(tau1-minlogtau)/
odpos1=min(real(NumTau,dp),max(0.0_dp,1.0+(tau1-minlogtau)/dlogtau))
odpos1=min(real(NumTau,dp),max(0.0_dp,1.0_dp+(tau1-minlogtau)/dlogtau))
iodpo1=int(odpos1)
dodpo1=odpos1-real(iodpo1,dp)
iodp11=min(NumTau,iodpo1+1)
Expand All @@ -548,7 +548,7 @@ subroutine photoion (phi,hcolum_in,hcolum_out,vol)
! find the table positions for the optical depth (outgoing)
tau1=log10(max(1.0e-20_dp,tauh_out))
! odpos1=min(1.0d0*NumTau,max(0.0d0,1.0d0+(tau1-minlogtau)/
odpos1=min(real(NumTau,dp),max(0.0_dp,1.0+(tau1-minlogtau)/dlogtau))
odpos1=min(real(NumTau,dp),max(0.0_dp,1.0_dp+(tau1-minlogtau)/dlogtau))
iodpo1=int(odpos1)
dodpo1=odpos1-real(iodpo1)
iodp11=min(NumTau,iodpo1+1)
Expand Down