ran1.F

Go to the documentation of this file.

#include "cppdefs.h"
      SUBROUTINE ran1_s (harvest)
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  Lagged  Fibonacci  generator  combined with two  Marsaglia  shift   !
!  sequences. On output, returns as HARVEST a uniform random deviate   !
!  between  0.0  and  1.0  (exclusive of the endpoint values).  This   !
!  generator has the same calling and  initialization conventions as   !
!  Fortran 90 RANDOM_NUMBEER routine.  Use RAN_SEED to initialize or   !
!  reinitialize a particular sequence.  The period of this generator   !
!  is about 8.5E+37, and it fully vectorizes.  Validy of the integer   !
!  model assumend by this generator is tested at initialization.       !
!                                                                      !
!  Scalar version adapted from Numerical recipes.                      !
!                                                                      !
!  Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery,    !
!     1996:  Numerical Recipes in Fortran 90,  The Art of Parallel     !
!     Scientific Computing, 2nd Edition, Cambridge Univ. Press.        !
!                                                                      !
!=======================================================================
!
      USE mod_kinds
      USE ran_state, ONLY : ran_init
      USE ran_state, ONLY : iran0, jran0, kran0, nran0, mran0
      USE ran_state, ONLY : amm, lenran, rans
!
!  Imported variable declarations.
!
      real(r8), intent(out) :: harvest
!
!-----------------------------------------------------------------------
!  Compute an uniform random deviate (scalar version).
!-----------------------------------------------------------------------
!
!  Initialize.
!
      IF (lenran.lt.1) CALL ran_init (1_i8b)
!
!  Update Fibonacci generator, which has a period p*p+p+1 (p=2^(31)-69).
!
      rans=iran0-kran0
      IF (rans.lt.0) rans=rans+2147483579_i8b
      iran0=jran0
      jran0=kran0
      kran0=rans
!
!  Update Marsaglia shift sequence.
!
      nran0=ieor(nran0,ishft(nran0,13))
      nran0=ieor(nran0,ishft(nran0,-17))
      nran0=ieor(nran0,ishft(nran0,5))
!
!  Once only per cycle, advance sequence by 1, shortening its period to
!  2^(32)-2.
!
      IF (nran0.eq.1) nran0=270369_i8b
!
!  Update Marsaglia shift sequence with perios 2^(32)-1.
!
      mran0=ieor(mran0,ishft(mran0,5))
      mran0=ieor(mran0,ishft(mran0,-13))
      mran0=ieor(mran0,ishft(mran0,6))
!
!  Wrap=around addition.
!
      rans=ieor(nran0,rans)+mran0
!
!  Make the results positive definite (note that AMM is negative).
!
      harvest=amm*merge(rans,not(rans), rans<0)
 
      RETURN
      END SUBROUTINE ran1_s
 
      SUBROUTINE ran1_v (harvest)
!
!=======================================================================
!                                                                      !
!  Lagged  Fibonacci  generator  combined with two  Marsaglia  shift   !
!  sequences. On output, returns as HARVEST a uniform random deviate   !
!  between  0.0  and  1.0  (exclusive of the endpoint values).  This   !
!  generator has the same calling and  initialization conventions as   !
!  Fortran 90 RANDOM_NUMBEER routine.  Use RAN_SEED to initialize or   !
!  reinitialize a particular sequence.  The period of this generator   !
!  is about 8.5E+37, and it fully vectorizes.  Validy of the integer   !
!  model assumend by this generator is tested at initialization.       !
!                                                                      !
!  Vector version adapted from Numerical recipes.                      !
!                                                                      !
!  Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery,    !
!     1996:  Numerical Recipes in Fortran 90,  The Art of Parallel     !
!     Scientific Computing, 2nd Edition, Cambridge Univ. Press.        !
!                                                                      !
!=======================================================================
!
      USE mod_kinds
      USE ran_state, ONLY : ran_init
      USE ran_state, ONLY : iran, jran, kran, nran, mran
      USE ran_state, ONLY : amm, lenran, ranv
!
!  Imported variable declarations.
!
      real(r8), dimension(:), intent(out) :: harvest
!
!  Local variable declarations.
!
      integer(i8b) :: n
!
!-----------------------------------------------------------------------
!  Compute an uniform random deviate (scalar version).
!-----------------------------------------------------------------------
!
!  Initialize.
!
      n=SIZE(harvest)
      IF (lenran.lt.n+1) CALL ran_init (n+1_i8b)
!
!  Update Fibonacci generator, which has a period p*p+p+1 (p=2^(31)-69).
!
      ranv(1:n)=iran(1:n)-kran(1:n)
      WHERE (ranv(1:n).lt.0)                                            &
     &  ranv(1:n)=ranv(1:n)+2147483579_i8b
      iran(1:n)=jran(1:n)
      jran(1:n)=kran(1:n)
      kran(1:n)=ranv(1:n)
!
!  Update Marsaglia shift sequence.
!
      nran(1:n)=ieor(nran(1:n),ishft(nran(1:n),13))
      nran(1:n)=ieor(nran(1:n),ishft(nran(1:n),-17))
      nran(1:n)=ieor(nran(1:n),ishft(nran(1:n),5))
!
!  Once only per cycle, advance sequence by 1, shortening its period to
!  2^(32)-2.
!
      WHERE (nran(1:n).eq.1)                                            &
     &  nran(1:n)=270369_i8b
!
!  Update Marsaglia shift sequence with perios 2^(32)-1.
!
      mran(1:n)=ieor(mran(1:n),ishft(mran(1:n),5))
      mran(1:n)=ieor(mran(1:n),ishft(mran(1:n),-13))
      mran(1:n)=ieor(mran(1:n),ishft(mran(1:n),6))
!
!  Wrap=around addition.
!
      ranv(1:n)=ieor(nran(1:n),ranv(1:n))+mran(1:n)
!
!  Make the results positive definite (note that AMM is negative).
!
      harvest=amm*merge(ranv(1:n),not(ranv(1:n)), ranv(1:n)<0 )
 
      RETURN
      END SUBROUTINE ran1_v