Optimizing a Fortran subroutine

I've written a minimal implementation for the fast xoroshiro128plus pseudo-random number generator in Fortran to replace the intrinsic random_number. This implementation is quite fast (4X faster than random_number) and the quality is good enough for my purposes, I don't use it for cryptography applications.

My question is how can I optimize this subroutine to get the last drop of performance from my compiler, even 10% improvement is appreciated. This subroutine is to be used in tight loops inside long simulations. I'm interested more in generating a single random number at a time and not big vectors or nD arrays at once.

Here is a test program to give you some context:

program test_xoroshiro128plus
   implicit none
   integer, parameter :: n = 10000
   real*8  :: A(n,n)
   integer :: i, j, t0, t1, count_rate, count_max

   call system_clock(t0, count_rate, count_max)
   do j = 1,n
      do i = 1,n
         call drand128(A(i,j))
      end do
   end do
   ! call drand128(A)  # works also with 2D 
   call system_clock(t1)

   print *, "Time :", real(t1-t0)/count_rate
   print *, "Mean :", sum(A)/size(A), char(10), A(1:2,1:3)

 contains

   impure elemental subroutine drand128(r)
      real*8, intent(out) :: r
      integer*8 :: s0 = 113, s1 = 19937
      s1 = xor(s0,s1)
      s0 = xor(xor(ior(ishft(s0,55), ishft(s0,-9)),s1), ishft(s1,14))
      s1 = ior(ishft(s1,36), ishft(s1,-28))
      r = ishft(s0+s1, -1) / 9223372036854775808.d0
   end 

end program