Thank you very much! @Arjen

It is loop over rnor() is not bad. Actually the rnor() seems faster than many other RNGs. I am just curious if there are faster ways to directly generator n gaussian random numbers instead of looping over rnor() for n times.

I have a subroutine

```
subroutine rk4_ti_fullvec_test ( x0, np, nstep, q, h, fi_gi_in, nd, x )
! https://stackoverflow.com/questions/69147944/is-there-room-to-further-optimize-the-stochastic-rk-fortran-90-code
! https://stackoverflow.com/questions/32809769/how-to-pass-subroutine-names-as-arguments-in-fortran
use random
use fg
implicit none
integer(kind = i8), intent(in) :: np
integer(kind = i8), intent(in) :: nstep
integer(kind = i8), intent(in) :: nd
procedure(fi_gi_fullvec_03) :: fi_gi_in
real(kind = r8), intent(in) :: q,h
real(kind = r8), intent(in) :: x0(np,nd)
real ( kind = r8 ), parameter :: alphas(4) = [ 0.47012396888046_r8 &
,0.36597075368373_r8 &
,0.08906615686702_r8 &
,0.07483912056879_r8 ]
real(kind = r8), parameter :: as(4,4) = reshape([ &
& 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, &
& 2.71644396264860_r8, 0.0_r8, 0.0_r8, 0.0_r8, &
& -6.95653259006152_r8, 0.78313689457981_r8, 0.0_r8, 0.0_r8, &
& 0.0_r8, 0.48257353309214_r8, 0.26171080165848_r8, 0.0_r8 ], [4,4])
real(kind = r8), parameter :: qs(4) = [ 2.12709852335625_r8, &
& 2.73245878238737_r8, &
& 11.22760917474960_r8, &
& 13.36199560336697_r8 ]
real(kind = r8) :: ks(np,nd,4)
real(kind = r8) :: xs(np,nd,4)
integer(kind = i8) :: i,j,k,l,m,n
real(kind = r8), intent(out) :: x(np,nd,0:nstep)
real(kind = r8) :: xstar(np,nd)
real ( kind = r8 ) :: f(np,nd), g(np,nd)
real(kind = r8) :: warray(np,nd,4)
real(kind = r8) :: sigma(4)
!real (kind = r8), allocatable :: normal(:,:,:,:)
!allocate(normal(np,nd,4,nstep))
!normal = reshape(rnor_vec(int(np*nd*4*nstep,kind=i4)),shape(normal))
sigma=sqrt(qs*q/h)
x(:,:,0) = x0
do k = 1, nstep
xstar = x(:,:,k-1)
do j = 1,4
do l = 1,nd
xs(:,l,j) = x(:,l,k-1) + matmul(ks(:,l,:j-1),as(:j-1,j))
do n = 1,np
warray(n,l,j) = rnor()*sigma(j) !gaussian(2)*sqrt(qs(j)*qoh)
enddo
enddo
call fi_gi_in( np, nd, xs(:,:,j), x0(1,:), f, g )
ks(:,:,j) = h * ( f + g*warray(:,:,j) )
!ks(:,:,j) = h * ( f + g*normal(:,:,j,k)*sigma(j) )
xstar = xstar + alphas(j)*ks(:,:,j)
enddo
x(:,:,k) = xstar
enddo
return
end subroutine rk4_ti_fullvec_test
pure subroutine fi_gi_fullvec_03 (np,nd, x, x0, f, g )
use constants
implicit none
integer (kind = i8), intent(in) :: np,nd
real ( kind = r8 ), intent(in) :: x(np,nd)
real ( kind = r8 ), intent(in) :: x0(nd)
real ( kind = r8 ), intent(out) :: f(np,nd)
real ( kind = r8 ), intent(out) :: g(np,nd)
integer (kind = i8) :: i
real ( kind = r8 ), parameter :: gp(2) = [one,zero]
f(:,1) = -x(:,2)*x(:,1)
f(:,2) = -x(:,2) + x0(2)
g(:,1) = gp(1)
g(:,2) = gp(2)
return
end subroutine fi_gi_fullvec_03
```

Inside the loop of subroutine rk4_ti_fullvec_test, I need to generator many gaussian random number such as

```
do n = 1,np
warray(n,l,j) = rnor()*sigma(j) !gaussian(2)*sqrt(qs(j)*qoh)
enddo
```

where np can be 10^6, nstep can be 1000, nd=2, j=1:4. So the speed of rnor() becomes a little bit crucial since need to generate like 10^6*1000*4*2 random gaussian numbers.

If there are ways to directly generator many gaussian random number and store them in an array outside the loops, this may be efficient.