Skip to content

Commit

Permalink
Use double precision for noise generator.
Browse files Browse the repository at this point in the history
  • Loading branch information
@SamuelTrahanNOAA
SamuelTrahanNOAA committed May 6, 2022
1 parent 290fb0f commit 471baa8
Show file tree
Hide file tree
Showing 2 changed files with 23 additions and 23 deletions.
44 changes: 22 additions & 22 deletions mersenne_twister.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -160,7 +160,7 @@
!
!$$$
module mersenne_twister
use kinddef, only: kind_phys
use kinddef, only: kind_dbl_prec
private
! Public declarations
public random_stat
Expand Down Expand Up @@ -189,7 +189,7 @@ module mersenne_twister
integer:: mti=n+1
integer:: mt(0:n-1)
integer:: iset
real(kind_phys):: gset
real(kind_dbl_prec):: gset
end type
! Saved data
type(random_stat),save:: sstat
Expand Down Expand Up @@ -297,8 +297,8 @@ subroutine random_setseed_t(inseed,stat)
!> This function generates random numbers in functional mode.
function random_number_f() result(harvest)
implicit none
real(kind_phys):: harvest
real(kind_phys) :: h(1)
real(kind_dbl_prec):: harvest
real(kind_dbl_prec) :: h(1)
if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
call random_number_t(h,sstat)
harvest=h(1)
Expand All @@ -307,7 +307,7 @@ function random_number_f() result(harvest)
!> This subroutine generates random numbers in interactive mode.
subroutine random_number_i(harvest,inseed)
implicit none
real(kind_phys),intent(out):: harvest(:)
real(kind_dbl_prec),intent(out):: harvest(:)
integer,intent(in):: inseed
type(random_stat) stat
call random_setseed_t(inseed,stat)
Expand All @@ -317,15 +317,15 @@ subroutine random_number_i(harvest,inseed)
!> This subroutine generates random numbers in saved mode; overloads Fortran 90 standard.
subroutine random_number_s(harvest)
implicit none
real(kind_phys),intent(out):: harvest(:)
real(kind_dbl_prec),intent(out):: harvest(:)
if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
call random_number_t(harvest,sstat)
end subroutine
! Subprogram random_number_t
!> This subroutine generates random numbers in thread-safe mode.
subroutine random_number_t(harvest,stat)
implicit none
real(kind_phys),intent(out):: harvest(:)
real(kind_dbl_prec),intent(out):: harvest(:)
type(random_stat),intent(inout):: stat
integer j,kk,y
integer tshftu,tshfts,tshftt,tshftl
Expand Down Expand Up @@ -353,9 +353,9 @@ subroutine random_number_t(harvest,stat)
y=ieor(y,iand(tshftt(y),tmaskc))
y=ieor(y,tshftl(y))
if(y.lt.0) then
harvest(j)=(real(y,kind_phys)+2.0_kind_phys**32)/(2.0_kind_phys**32-1.0_kind_phys)
harvest(j)=(real(y,kind_dbl_prec)+2.0_kind_dbl_prec**32)/(2.0_kind_dbl_prec**32-1.0_kind_dbl_prec)
else
harvest(j)=real(y,kind_phys)/(2.0_kind_phys**32-1.0_kind_phys)
harvest(j)=real(y,kind_dbl_prec)/(2.0_kind_dbl_prec**32-1.0_kind_dbl_prec)
endif
stat%mti=stat%mti+1
enddo
Expand All @@ -364,8 +364,8 @@ subroutine random_number_t(harvest,stat)
!> This subrouitne generates Gaussian random numbers in functional mode.
function random_gauss_f() result(harvest)
implicit none
real(kind_phys):: harvest
real(kind_phys) :: h(1)
real(kind_dbl_prec):: harvest
real(kind_dbl_prec) :: h(1)
if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
call random_gauss_t(h,sstat)
harvest=h(1)
Expand All @@ -374,7 +374,7 @@ function random_gauss_f() result(harvest)
!> This subrouitne generates Gaussian random numbers in interactive mode.
subroutine random_gauss_i(harvest,inseed)
implicit none
real(kind_phys),intent(out):: harvest(:)
real(kind_dbl_prec),intent(out):: harvest(:)
integer,intent(in):: inseed
type(random_stat) stat
call random_setseed_t(inseed,stat)
Expand All @@ -384,18 +384,18 @@ subroutine random_gauss_i(harvest,inseed)
!> This subroutine generates Gaussian random numbers in saved mode.
subroutine random_gauss_s(harvest)
implicit none
real(kind_phys),intent(out):: harvest(:)
real(kind_dbl_prec),intent(out):: harvest(:)
if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
call random_gauss_t(harvest,sstat)
end subroutine
! Subprogram random_gauss_t
!> This subroutine generates Gaussian random numbers in thread-safe mode.
subroutine random_gauss_t(harvest,stat)
implicit none
real(kind_phys),intent(out):: harvest(:)
real(kind_dbl_prec),intent(out):: harvest(:)
type(random_stat),intent(inout):: stat
integer mx,my,mz,j
real(kind_phys) :: r2(2),r,g1,g2
real(kind_dbl_prec) :: r2(2),r,g1,g2
mz=size(harvest)
if(mz.le.0) return
mx=0
Expand Down Expand Up @@ -430,14 +430,14 @@ subroutine random_gauss_t(harvest,stat)
contains
!> This subroutine contains numerical Recipes algorithm to generate Gaussian random numbers.
subroutine rgauss(r1,r2,r,g1,g2)
real(kind_phys),intent(in):: r1,r2
real(kind_phys),intent(out):: r,g1,g2
real(kind_phys) :: v1,v2,fac
v1=2.*r1-1._kind_phys
v2=2.*r2-1._kind_phys
real(kind_dbl_prec),intent(in):: r1,r2
real(kind_dbl_prec),intent(out):: r,g1,g2
real(kind_dbl_prec) :: v1,v2,fac
v1=2.*r1-1._kind_dbl_prec
v2=2.*r2-1._kind_dbl_prec
r=v1**2+v2**2
if(r.lt.1.) then
fac=sqrt(-2._kind_phys*log(r)/r)
fac=sqrt(-2._kind_dbl_prec*log(r)/r)
g1=v1*fac
g2=v2*fac
endif
Expand Down Expand Up @@ -483,7 +483,7 @@ subroutine random_index_t(imax,iharvest,stat)
type(random_stat),intent(inout):: stat
integer,parameter:: mh=n
integer i1,i2,mz
real(kind_phys) :: h(mh)
real(kind_dbl_prec) :: h(mh)
mz=size(iharvest)
do i1=1,mz,mh
i2=min((i1-1)+mh,mz)
Expand Down
2 changes: 1 addition & 1 deletion stochy_patterngenerator.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -243,7 +243,7 @@ subroutine getnoise(rpattern,noise_e,noise_o)
real(kind_phys), intent(out) :: noise_o(len_trio_ls,2)
! generate white noise with unit variance in spectral space
type(random_pattern), intent(inout) :: rpattern
real :: noise(2*ndimspec)
real(kind_dbl_prec) :: noise(2*ndimspec)
integer nm,nn
call random_gauss(noise,rpattern%rstate)
noise(1) = 0.; noise(ndimspec+1) = 0.
Expand Down

0 comments on commit 471baa8

Please sign in to comment.