!=======================================================================
!
!  sderk90 module
!
!  Version 2.1
!
!  Main interface subroutine:
!
!     sderk (t, tout)
!                      
!  This subroutine generates strong solutions to Ito stochastic 
!  differential equations of the form
!    
!     dy = a(y, t) dt + b(y, t) dW
!
!  where y is a vector, a and b are vector-valued functions, and W is
!  a scalar Wiener process.  This routine implements Runge-Kutta
!  methods of order 1 and 1.5.
!  This routine also constructs a Brownian tree to realize the same
!  Brownian path, independent of the order of the method or the 
!  internal stepsize, as long as the routine is called with the same
!  output times.
!          
!  Random numbers are generated using an internal, portable generator,
!  so calling other generators between integration steps will not
!  interfere with the integrator.
!
!
!  References:
!    
!  Pamela Marion Burrage, "Runge-Kutta Methods for Stochastic
!    Differential Equations," Ph.D. thesis, University of Queensland
!    (1999).
!    
!  P. E. Kloeden and E. Platen, "Higher-Order Implicit Strong Numerical
!    Schemes for Stochastic Differential Equations," J. Stat. Phys. 66,
!    283 (1992).
!    
!  Peter E. Kloeden and Eckhard Platen, _Numerical Solution of
!    Stochastic Differential Equations_ (Springer, Berlin, 1992).
!    
!  J. G. Gaines and T. J. Lyons, "Variable Step Size Control in the
!    Numerical Solution of Stochastic Differential Equations,"
!    SIAM J. Appl. Math. 57, 1455 (1997).
!
!                      
!  Interface:
!
!  This module defines:
!
!    sderk_int_type -> (integer parameter) the integer type used
!                  by the internal random-number generator; this
!                  parameter is needed to declare some integer
!                  scalars and arrays of type integer(sderk_int_type)
!                  to inteface with the internal rng (e.g., to seed
!                  it or to save/restore its state).
!
!  Most of the communication with this integrator is through a
!    user-defined module called "sderk_support". A sample module
!    is included as a separate file with this software.  This module
!    should define several parameters, variables, and functions:
!                    
!    sderk_prec -> (default integer parameter) the real/complex kind
!                  for the precision type to use for all computations,
!                  e.g., this should be the output of selected_real_kind.
!
!    sderk_mf  -> (default integer) method flag to set which integration
!                 method to use (see below for available options).
!                      
!    sderk_rand_mf -> (default integer) method flag to set which
!                 method for generating random numbers to use.  The
!                 allowed values are documented in the random_pl
!                 (version 2.0) package by the same author, but 
!                 a reasonable choice is 301 (or 204, to switch
!                 methods to see if the random-number generator is
!                 causing funny results).
!
!    sderk_tol -> (real(sderk_prec)) tolerance for iterative solution,
!                 this is needed for implicit methods and ignored
!                 otherwise (but still needs to be defined). A  
!                 sensible number is 1e-10 for double precision.
!                    
!    sderk_y -> (real/complex array) the solution array, where
!               the type/dimension are set according to the
!               user's requirements.  Because the initial
!               value is set by the user, the user must control
!               the allocation for this array.  May be declared
!               as a pointer to a more sensible array name.
!    
!    sderk_a      ->  \
!    sderk_b      ->   |  Always allocate these
!    sderk_c      ->   |
!    sderk_ybar   ->  /
!    sderk_aybarp ->  \
!    sderk_aybarm ->   |
!    sderk_bybarp ->   |
!    sderk_bybarm ->   | For any order 1.5 stochastic method
!    sderk_aphip  ->   | Also any order 4 deterministic (phi's only)
!    sderk_aphim  ->   |  
!    sderk_bphip  ->   |
!    sderk_bphim  ->  /
!                   auxilliary storage arrays, of the same type and
!                   dimensions as sderk_solution.  These should be   
!                   declared and allocated by the user.  It is easiest
!                   to just define and allocate all of these, but in
!                   some cases, as marked above, not all of these are
!                   actually needed; the unneeded ones can be declared
!                   "allocatable" and not actually allocated, if the
!                   user is memory-starved.
!    sderk_func_a -> subroutine of the form
!                    subroutine sderk_func_a(t, y, ydot)
!                      real(sderk_prec), intent(in) :: t
!                      <solution type>, intent(in)  :: y
!                      <solution type>, intent(out) :: ydot
!                    which returns the derivative array a(y, t). Here,
!                    <solution type> is the array type/dimension
!                    defined above for the solution arrays, e.g., this
!                    could be "real(sderk_prec), dimension(100,3)".
!
!    sderk_func_b -> subroutine of the form
!                    subroutine sderk_func_b(t, y, ydot)
!                      real(sderk_prec), intent(in) :: t
!                      <solution type>, intent(in)  :: y
!                      <solution type>, intent(out) :: ydot
!                    which returns the derivative array b(y, t).
!
!
!  Calling interface:
!
!  The integrator is called to advance in time via
!      call sderk (t, tout)
!    where
!
!    t -> (real(sderk_prec)) initial value of the time parameter;
!       on successful integration, this is updated to tout
!
!    tout -> (real(sderk_prec)) the desired output time
!
!
!  Before calling this routine, the integrator must be initialized by
!    a call of the form
!      call init_sderk(nsteps, seed)
!    where
!
!    nsteps -> (default integer) number of steps to take from t to tout;
!            must be a power of 2 (1, 2, 4, 8, etc.).  Two runs with
!            the same t and tout but different values of nsteps will
!            be realizations of the same Brownian path.
!
!    seed -> (optional integer(sderk_int_type)) seed to initialize
!            random number generators; a default seed is used if
!            this seed is omitted.
!
!    Note that if this routine is called again with other values of
!    nsteps (i.e., nsteps is not constant over a long integration),
!    the results of several runs with different nsteps values
!    may not be on the same Brownian path.  The only way to
!    ensure that a consistent Brownian path is used is to keep
!    nsteps fixed over a given run.  In the future, an adaptive
!    driver will use this facility to modify nsteps, but saving
!    the integrator state and resynchronizing the various
!    random-number generators to enforce consistent Brownian paths.
!
!
!  Numerical methods:
!  
!      This package supports a number of integration methods, set by
!        the method flag mf, a two-digit integer.
!        The first digit is a flag (1 or 2), with 1 indicating
!        explicit and 2 implicit evolution of the deterministic
!        part of the equation (i.e. a(y, t)).
!        Solution of the implicit equations is accomplished by
!        simple functional iteration; setting this digit to
!        3 instead of two causes the integrator to report the
!        number of iterations needed to converge.
!        The second digit is the order of the Runge-Kutta solver
!        Legal values are
!          11 (first-order explicit deterministic,
!                first-order stochastic)
!          12 (fourth-order explicit deterministic,
!                first-order stochastic)
!          13 (explicit order 1.5 RK method in Kloeden/Platen)
!          21 (second-order implicit deterministic,
!                first-order stochastic)
!          22 (fourth-order implicit deterministic,
!                first-order stochastic)
!          23 (implicit order 1.5 RK method in Kloeden/Platen; order
!                2 deterministic step)
!          24 (modification of 23 to have 4th-order implicit
!                deterministic step, but has possibly worse
!                than order 1.5 convergence)
!          25 (modification of 23 to have 6th-order implicit
!                deterministic step, also with possibly worse
!                than order 1.5 convergence)
! 
!                
!  Random numbers:
!          
!  This package uses an embedded, internal random-number generator, 
!    which should not be disturbed during normal operation.  For
!    checkpointing and restarting between calls, three subroutines
!    for accessing the generator states are available:
!
!    subroutine sderk_get_rand_state_size(dim1, dim2) ->
!             gets dimensions of the istate array needed for the
!             call to sderk_get_rand_state
!
!    subroutine sderk_get_rand_state(istate) ->
!             gets the state of the rng, putting the result in istate;
!             istate is an integer(sderk_int_type), allocatable array,
!             which is allocated to the proper size automatically.
!
!    subroutine sderk_set_rand_state(istate) ->
!             initializes and resets the generator to the state when
!             istate was returned by a call to sderk_get_rand_state.
!
!    Thus, the procedure to get the state of the rng would be:
!      integer :: dim1, dim2
!      integer(sderk_int_kind), dimension(:,:), allocatable :: istate
!      call sderk_get_rand_state_size(dim1, dim2)
!      allocate( istate(dim1, dim2) )
!      call sderk_get_rand_state(istate)
!
!    and a later call sderk_set_rand_state(istate) would restore
!    the generator to that state.
!
!
!  Also, there are two interface routines for retrieving the random 
!    numbers used for the last call to sderk:
!    
!    sderk_get_I1() -> function returning a vector of real(sderk_prec)
!                      numbers and of length nsteps of the Wiener  
!                      increments I1 (i.e., dW) used for the internal
!                      integration steps.
!    
!    sderk_get_I10() -> similar function, but returning the I10
!                       integrals instead.
!
!=======================================================================


module sderk90

  use sderk_support, only :                                            &
           sderk_func_a, sderk_func_b,                                 &
           sderk_prec, sderk_tol, sderk_mf, sderk_y,                   &
           sderk_a, sderk_b, sderk_c, sderk_ybar,                      &
           sderk_aybarp, sderk_aybarm, sderk_bybarp, sderk_bybarm,     &
           sderk_aphip, sderk_aphim, sderk_bphip, sderk_bphim,         &
           sderk_rand_mf


  private
  public :: sderk_int_kind, sderk_get_rand_state_size
  public :: sderk_get_rand_state, sderk_set_rand_state
  public :: init_sderk, sderk, sderk_get_I1, sderk_get_I10


  !!! local storage for integrator stuff

  ! aliases
  integer, parameter  :: sderk_int_kind = selected_int_kind(9)
  integer, parameter  :: rpk = sderk_int_kind     
  integer, parameter  :: wp = sderk_prec
  integer, parameter  :: mf = sderk_mf

  ! maximum number of iterations before convergence fails    
  integer, parameter :: maxit = 50
  real(wp)           :: tol

  ! storage for random deviates
  real(wp), dimension(:),   allocatable :: I1, I10
  real(wp), dimension(:,:), allocatable :: I1raw
  real(wp), dimension(2) :: z

  ! Brownian tree stuff
  integer :: nsteps, nsteps_last = -1, nsteps_most = - 1, treedepth

  !!! length of Brownian tree base, also the number of random numbers
  !!!   used to calculate the Ito integral I_10.
  !!!   Must be at least 2.
  integer, parameter :: tree_base = 20

  !!! select Brownian tree refinement method
  !!!   choices are levyarea for Levy area method (just refine I_1's)
  !!!   or splitboth for transformation method by Mauthner/Burrage/Burrage
  integer, parameter :: levyarea = 1, splitboth = 2
  integer, parameter :: refinement_method = splitboth

  !!! local storage for embedded random-number generator

  ! default seeds
  integer(rpk), parameter :: dseed1 = 310952
  integer(rpk), parameter :: dseed2 = 314159
  integer(rpk), parameter :: dseed3 = 4357

  ! parse out digits of sderk_rand_mf
  integer, parameter :: rand_pl_mf = sderk_rand_mf
  integer, parameter :: rand_pl_mf_pri = sderk_rand_mf/100
  integer, parameter :: rand_pl_mf_mix = sderk_rand_mf/10 - rand_pl_mf_pri*10
  integer, parameter :: rand_pl_mf_scr = sderk_rand_mf - (sderk_rand_mf/10)*10
  
  ! temporary storage for fibonacci generator
  integer(rpk), dimension(1009) :: aa
  
  ! define storage size parameters
  integer, parameter :: scramble_stor_sz = 64
  integer, parameter :: fib_stor_sz = 100 
  integer, parameter :: lecuyer_stor_sz = 6
  integer, parameter :: mersenne_stor_sz = 624
  integer, parameter :: rand_pl_state_sz = &
                scramble_stor_sz+fib_stor_sz+lecuyer_stor_sz+mersenne_stor_sz+5
  
  ! define storage for multiple generators
  integer(rpk), allocatable, dimension(:,:), target :: lecuyer_stor_m
  integer(rpk), allocatable, dimension(:,:), target :: scramble_stor_m
  integer(rpk), allocatable, dimension(:,:), target :: fib_stor_m
  integer(rpk), allocatable, dimension(:,:), target :: mersenne_stor_m
  integer(rpk), allocatable, dimension(:), target :: initialized_m
  integer(rpk), allocatable, dimension(:), target :: fib_ctr_m, bd_lastout_m, &
                                                     mersenne_ptr_m

  ! define storage pointers
  integer(rpk), dimension(:), pointer :: lecuyer_stor
  integer(rpk), dimension(:), pointer :: scramble_stor
  integer(rpk), dimension(:), pointer :: fib_stor
  integer(rpk), dimension(:), pointer :: mersenne_stor
  integer(rpk), pointer :: initialized
  integer(rpk), pointer :: fib_ctr, bd_lastout, mersenne_ptr

  ! to keep track of where pointers are associated
  integer :: curr_ptr = -1


contains


!-----------------------------------------------------------------------
!
!  main subroutine: sderk
!
!-----------------------------------------------------------------------

subroutine sderk(t, tout)

  implicit none

  real(wp), intent(inout) :: t
  real(wp), intent(in)    :: tout

  integer      :: j, k, p, q
  real(wp)     :: bigstep, h
  real(wp), dimension(nsteps) :: A10, sumarr

  real(wp), dimension(tree_base) :: tmp1, tmp2
  real(wp) :: btmp1, btmp2, btmp3, btmp4


  !!! Make sure initializations already happened
  if ( .not. allocated( initialized_m ) ) then
    write(0,*) 'Error (SDERK): integrator called before initialized'
    stop
  end if

  if ( refinement_method .eq. levyarea ) then

  !!! Precompute time steps
  bigstep = tout - t
  h = bigstep / real(tree_base, kind=wp)

  !!! Generate base layer of random numbers; store in part of rwork
  I1raw(1,:) = nrand_pl_vec_multi(1, tree_base) * sqrt(h)

  !!! Descend Brownian tree, only retaining finest layer     
  do j = 1, treedepth

    ! Rearrange values to create space for new deviates
    do k = 2**(j-1), 2, -1
      I1raw(2*k-1, :) = I1raw(k, :)
    end do

    ! get new deviates
    h = bigstep/real(tree_base*2**j, kind=wp)
    if ( j .eq. 1 ) then
      I1raw(2, :) = nrand_pl_vec_multi(2, tree_base)
    else
      do k = 1, tree_base
        I1raw(2:(2**j):2, k) = nrand_pl_vec_multi(1+j, 2**(j-1))
      end do
    end if

    ! refine I1 onto smaller grid
    do k = 1, 2**j, 2
      tmp1 = I1raw(k, :)   * 0.5_wp
      tmp2 = I1raw(k+1, :) * sqrt(0.5_wp * h)
      I1raw(k, :)   = tmp1 + tmp2
      I1raw(k+1, :) = tmp1 - tmp2
    end do

  end do

  !!! Calculate I1's and I10's via Levy areas A10
  h = bigstep/nsteps
  I1raw = transpose(reshape(I1raw, (/ tree_base, nsteps /)))
  I1 = sum(I1raw, dim=2)

  sumarr = 0.0_wp
  A10 = 0.0_wp
  do p = 2, tree_base
    sumarr = sumarr + I1raw(:, p-1)
    A10 = A10 + sumarr
  end do
  do q = 2, tree_base
    A10 = A10 - I1raw(:, q) * real(q-1, kind=wp)
  end do
  A10 = A10 * h/real(tree_base, kind=wp) * 0.5_wp
  I10 = A10 + 0.5_wp * I1 * h

  else if ( refinement_method .eq. splitboth ) then

  !!! Precompute time steps
  bigstep = tout - t
  h = bigstep
                
  !!! Generate base layer of random numbers; store in part of rwork
  z = nrand_pl_vec_multi(1, 2)
  I1(1) = z(1) * sqrt(h)
  I10(1) = h*sqrt(h)/2 * (z(1) + z(2)/sqrt(3.0_wp))
  
  !!! Descend Brownian tree, only retaining finest layer     
  do j = 1, treedepth
  
    ! Rearrange values to create space for new deviates
    do k = 2**(j-1), 2, -1
      I1(2*k-1) = I1(k)
      I10(2*k-1) = I10(k)
    end do
  
    ! refine I1 onto smaller grid
    h = bigstep / 2**(j-1) ! h is coarser step here
    do k = 1, 2**j, 2
      z = nrand_pl_vec_multi(1+j, 2)
      btmp1 = -sqrt(h)*z(2)/4 - I1(k)/4 + 3*I10(k)/(2*h)
      btmp2 = h*sqrt(h)*z(1)/(8*sqrt(3.0_wp)) - h*I1(k)/8 + I10(k)/2
      btmp3 = sqrt(h)*z(2)/4 + 5*I1(k)/4 - 3*I10(k)/(2*h)
      btmp4 = -h*sqrt(h)*z(1)/(8*sqrt(3.0_wp)) &
               + h*sqrt(h)*z(2)/8 + h*I1(k)/4 -I10(k)/4
      I1(k) = btmp1
      I1(k+1) = btmp3
      I10(k) = btmp2
      I10(k+1) = btmp4
    end do

  end do

  h = bigstep/nsteps

  end if


  !!! Set iteration tolerance, if needed
  if ( mf/10 .eq. 2 .or. mf/10 .eq. 3 ) then
    if ( sderk_tol .lt. 100_wp * epsilon(1.0_wp) ) then
      tol = 100_wp * epsilon(1.0_wp)
      write(0,*) 'Warning (SDERK): resetting tolerance to ', tol
    else
      tol = sderk_tol
    end if
  end if


  !!! Now do integration by appropriate method
  if ( mf .eq. 11 .or. mf .eq. 12 ) then
    do k = 1, nsteps
      call rk1_ito_onestep(t, h, I1(k))
      t = t + h
    end do
  else if ( mf .eq. 13 ) then
    do k = 1, nsteps
      call rk15_ito_onestep(t, h, I1(k), I10(k))
      t = t + h
    end do
  else if ( mf .eq. 21 .or. mf .eq. 31 .or. mf .eq. 22 .or. mf .eq. 32 ) then
    do k = 1, nsteps
      call rk1_ito_imp_onestep(t, h, I1(k))
      t = t + h
    end do
  else if ( mf .eq. 23 .or. mf .eq. 33 .or. mf .eq. 24 .or. mf .eq. 34 .or. &
            mf .eq. 25 .or. mf .eq. 35 ) then
    do k = 1, nsteps
      call rk15_ito_imp_onestep(t, h, I1(k), I10(k))
      t = t + h
    end do
  else
    write(0,*) 'Error (SDERK): invalid mf'
    stop
  end if

  return

end subroutine sderk


!-----------------------------------------------------------------------
!
!  subroutine init_sderk(set_nsteps, optional seed)
!
!  Initializes random-number generators for sderk integrator, 
!  so that the integrator will use a number of substeps equal to
!  set_nsteps per call to sderk.  Calling this function again during
!  an evolution will reset the random number generators if
!  nsteps is larger than before.  The optional argument seed
!  is used to initialize the random number generators, with distinct
!  seeds producing distinct random sequences.  A default seed is 
!  used if this argument is omitted.
!
!  The integer argument set_nsteps must be a power of 2.
!  The integer(sderk_int_kind) argument seed must be in the range
!    [0, 2**30 - 3].
!
!-----------------------------------------------------------------------

subroutine init_sderk(set_nsteps, seed)

  implicit none

  integer, intent(in)                :: set_nsteps
  integer(rpk), optional, intent(in) :: seed

  nsteps = set_nsteps

  !!! Calculate depth of Brownian tree (treedepth = 0 means just base layer)
  treedepth = int( log(real(nsteps, kind=wp)) / log(2.0_wp) + 0.5_wp )
  if ( 2**treedepth .ne. nsteps ) then
    write(0,*) 'Error (INIT_SDERK): nsteps must be power of 2'
    stop
  end if

  !!! Manage storage for random deviates
  if ( allocated(I1) .and. ( nsteps .ne. nsteps_last ) ) then
    deallocate(I1, I10, I1raw)
  else if ( .not. allocated(I1) ) then
    allocate(I1(nsteps), I10(nsteps), I1raw(nsteps, tree_base))
    nsteps_last = nsteps
  end if

  !!! reallocate random number generators, if more are needed
  if ( nsteps .gt. nsteps_most ) then
    if ( present(seed) ) then
      call allocate_rand_pl_multi(treedepth+1, seed1=seed)
    else
      call allocate_rand_pl_multi(treedepth+1)
    end if
  end if

  nsteps_most = max(nsteps_most, nsteps)

end subroutine init_sderk


!-----------------------------------------------------------------------
!
!  subroutine sderk_get_rand_state_size(dim1, dim2)
!
!  Gets the dimensions of an array needed for a call to
!  sderk_get_rand_state.
!
!-----------------------------------------------------------------------

subroutine sderk_get_rand_state_size(dim1, dim2)
  implicit none
  integer, intent(out) :: dim1, dim2
  if ( .not. allocated( initialized_m ) ) then
    write(0,*) 'Error (SDERK_GET_RAND_STATE_SIZE): generators not yet ', &
               'initialized'
    stop
  end if
  dim1 = rand_pl_state_sz
  dim2 = size(initialized_m)
end subroutine sderk_get_rand_state_size


!-----------------------------------------------------------------------
!
!  subroutine sderk_get_rand_state(istate)
!
!  Gets the state of the random number generators, stored in the
!  allocatable integer(sderk_int_type) array istate.  This array
!  will be allocated by this subroutine.
!
!-----------------------------------------------------------------------

subroutine sderk_get_rand_state(istate)

  implicit none

  integer(rpk), dimension(:,:), intent(out) :: istate
  integer :: j, dim1, dim2

  call sderk_get_rand_state_size(dim1, dim2)
  if ( size(istate, 1) .ne. dim1 ) then
    write(0,*) 'Error (SDERK_GET_RAND_STATE): istate has wrong first dim'
    stop
  end if
  if ( size(istate, 2) .ne. dim2 ) then
    write(0,*) 'Error (SDERK_GET_RAND_STATE): istate has wrong second dim'
    stop
  end if

  do j = 1, size(istate, 2)
    call get_rand_pl_state_multi(j, istate(:,j))
  end do

end subroutine sderk_get_rand_state


!-----------------------------------------------------------------------
!
!  subroutine sderk_set_rand_state(istate)
!
!  Resets the state of the random number generators, stored in the
!  integer(sderk_int_type) array istate, which was retrieved from
!  a previous call to sderk_get_rand_state.  This subroutine
!  will also initialize the random number generators, so a prior
!  call to init_sderk is unnecessary.
!
!-----------------------------------------------------------------------

subroutine sderk_set_rand_state(istate)

  implicit none

  integer(rpk), dimension(:,:), intent(in) :: istate
  integer :: j

  if ( size(istate, 1) .ne. rand_pl_state_sz ) then
    write(0,*) 'Error (SDERK_SET_RAND_STATE): istate has wrong first dim'
    stop
  end if

  call init_sderk( 2**( size(istate, 2) - 1 ) )

  do j = 1, size(istate, 2)
    call set_rand_pl_state_multi(j, istate(:,j))
  end do

end subroutine sderk_set_rand_state


!-----------------------------------------------------------------------
!
!  subroutine rk1_ito_onestep(t, h, dW)
!
!  Calculates one order 1.0 stochastic Runge-Kutta iteration for an
!  Ito SDE. The parameter t is the time at the
!  beginning of the interval of length h.
!  dW is the Wiener increment during this interval.
!  This subroutine requires the user-defined arrays sderk_a, sderk_b,
!  and sderk_c.
!
!  The method here is that of Kloeden and Platen.
!
!-----------------------------------------------------------------------

subroutine rk1_ito_onestep(t, h, dW)

  implicit none     

  real(wp), intent(in) :: t, h, dW


  !!! Compute stochastic increment and store in sderk_c
  call sderk_func_a(t, sderk_y, sderk_a)
  call sderk_func_b(t, sderk_y, sderk_b)
  sderk_ybar = sderk_y + sderk_a * h + sderk_b * sqrt(h)

  call sderk_func_b(t, sderk_ybar, sderk_c)
  sderk_c = sderk_b * dW + (sderk_c - sderk_b) * (dW*dW-h)/(2.0_wp*sqrt(h))


  !!! Compute deterministic increment
  if ( mf .eq. 11 ) then
    sderk_y = sderk_y + sderk_a * h
  else if ( mf .eq. 12 ) then
    call rk4_onestep(t, h)
  end if

  sderk_y = sderk_y + sderk_c

  return

end subroutine rk1_ito_onestep


!-----------------------------------------------------------------------
!
!  subroutine rk1_ito_imp_onestep(t, h, dW)
!
!  Calculates one implicit order 1.0 stochastic Runge-Kutta iteration
!  for an Ito SDE. The parameter t is the time at the
!  beginning of the interval of length h.
!  dW is the Wiener increment during this interval.
!  This subroutine requires the user-defined storage arrays sderk_a,
!  sderk_b, sderk_c, sderk_ybar, sderk_aphip, sderk_aphim,
!  sderk_bphip, and sderk_bphim.
!
!  The method here is that of Kloeden and Platen, with degree of
!  implicitness 1/2.
!
!-----------------------------------------------------------------------

subroutine rk1_ito_imp_onestep(t, h, dW)

  implicit none

  real(wp), intent(in) :: t, h, dW


  !!! Compute stochastic increment and store in sderk_c
  call sderk_func_a(t, sderk_y, sderk_a)       
  call sderk_func_b(t, sderk_y, sderk_b)
  sderk_ybar = sderk_y + sderk_a * h + sderk_b * sqrt(h)

  call sderk_func_b(t, sderk_ybar, sderk_c)
  sderk_c = sderk_b * dW + (sderk_c - sderk_b) * (dW*dW-h)/(2.0_wp*sqrt(h))


  if ( (mf - 10*(mf/10)) .eq. 1 ) then
    sderk_c = sderk_c + 0.5_wp * sderk_a * h
    call ord_2_iterator(t, h)
  else
    call ord_4_iterator(t, h)
  end if

  return 

end subroutine rk1_ito_imp_onestep


!-----------------------------------------------------------------------
!
!  subroutine rk15_ito_onestep(t, h, dW, I10)
!
!  Calculates one order 1.5 stochastic Runge-Kutta iteration for an
!  Ito SDE. The parameter t is the time at the
!  beginning of the interval of length h.  dW and I10 are the
!  respective I1 and I10 integral values for this interval.
!  This routine requires the user-defined storage arrays
!  sderk_a, sderk_b, sderk_c, sderk_aybarp, sderk_aybarm,
!  sderk_bybarp, sderk_bybarm, sderk_aphip, sderk_aphim,
!  sderk_bphip, sderk_bphim, and sderk_ybar.
!
!  The method here is from Kloeden/Platen.
!
!-----------------------------------------------------------------------

subroutine rk15_ito_onestep(t, h, dW, I10)

  implicit none

  real(wp), intent(in) :: t, h, dW, I10


  ! compute a(y), b(y)
  call sderk_func_a(t, sderk_y, sderk_a)
  call sderk_func_b(t, sderk_y, sderk_b)

  ! compute the ybar(+) supporting value, and eval functions
  sderk_ybar = sderk_y + sderk_a * h + sderk_b * sqrt(h)
  call sderk_func_a(t + h, sderk_ybar, sderk_aybarp)
  call sderk_func_b(t + h, sderk_ybar, sderk_bybarp)

  ! compute the ybar(-) supporting value (store as c), and eval functions
  sderk_c = sderk_y + sderk_a * h - sderk_b * sqrt(h)
  call sderk_func_a(t - h, sderk_c, sderk_aybarm)
  call sderk_func_b(t - h, sderk_c, sderk_bybarm)

  ! compute the phi(+) supporting value (store as c), and eval functions
  sderk_c = sderk_ybar + sderk_bybarp * sqrt(h)
  call sderk_func_a(t + h, sderk_c, sderk_aphip)
  call sderk_func_b(t + h, sderk_c, sderk_bphip)

  ! compute the phi(-) supporting value (store as c), and eval functions
  sderk_c = sderk_ybar - sderk_bybarp * sqrt(h)
  call sderk_func_a(t + h, sderk_c, sderk_aphim)
  call sderk_func_b(t + h, sderk_c, sderk_bphim)

  ! compute the final result
  sderk_y = sderk_y + sderk_b * dW                                     &
           + (sderk_aybarp-sderk_aybarm)*I10/(2.0_wp*sqrt(h))          &
           + (sderk_aybarp+2.0_wp*sderk_a+sderk_aybarm)*h/4.0_wp       &
           + (sderk_bybarp-sderk_bybarm)*(dW*dW-h)/(4.0_wp*sqrt(h))    &
           + (sderk_bybarp-2.0_wp*sderk_b+sderk_bybarm)                &
             * (dW*h-I10)/(2.0_wp*h)                                   &
           + (sderk_bphip-sderk_bphim-sderk_bybarp+sderk_bybarm)       &
             * (dW*dW/3.0_wp-h)*dW/(4.0_wp*h)

  return

end subroutine rk15_ito_onestep


!-----------------------------------------------------------------------
!
!  subroutine rk15_ito_imp_onestep(t, h, dW, I10)
!
!  Calculates one (implicit) order 1.5 stochastic Runge-Kutta iteration
!  for an Ito SDE. The parameter t is the time at the
!  beginning of the interval of length h.  dW and I10 are the
!  respective I1 and I10 integral values for this interval.
!  This routine requires the user-defined arrays
!  sderk_a, sderk_b, sderk_c, sderk_aybarp, sderk_aybarm,
!  sderk_bybarp, sderk_bybarm, sderk_aphip, sderk_aphim,   
!  sderk_bphip, sderk_bphim, sderk_ybar.
!
!  The method here is from Kloeden/Platen, with degree of implicitness    
!  1/2.
!
!-----------------------------------------------------------------------

subroutine rk15_ito_imp_onestep(t, h, dW, I10)

  implicit none

  real(wp), intent(in) :: t, h, dW, I10


  call sderk_func_a(t, sderk_y, sderk_a)
  call sderk_func_b(t, sderk_y, sderk_b)

  ! compute the ybar(+) supporting value, and eval functions
  sderk_ybar = sderk_y + sderk_a * h + sderk_b * sqrt(h)
  call sderk_func_a(t + h, sderk_ybar, sderk_aybarp)
  call sderk_func_b(t + h, sderk_ybar, sderk_bybarp)

  ! compute the ybar(-) supporting value (store as c), and eval functions
  sderk_c = sderk_y + sderk_a * h - sderk_b * sqrt(h)
  call sderk_func_a(t - h, sderk_c, sderk_aybarm)
  call sderk_func_b(t - h, sderk_c, sderk_bybarm)

  ! compute the phi(+) supporting value (store as c), and eval functions
  sderk_c = sderk_ybar + sderk_bybarp * sqrt(h)
  call sderk_func_a(t + h, sderk_c, sderk_aphip)         
  call sderk_func_b(t + h, sderk_c, sderk_bphip)

  ! compute the phi(-) supporting value (store as c), and eval functions
  sderk_c = sderk_ybar - sderk_bybarp * sqrt(h)
  call sderk_func_a(t + h, sderk_c, sderk_aphim)
  call sderk_func_b(t + h, sderk_c, sderk_bphim)

  ! compute everything but the implicit part, store in c
  sderk_c = sderk_b * dW                                             &
           + (sderk_aybarp-sderk_aybarm)*(I10-0.5_wp*dW*h)           &
             / (2.0_wp*sqrt(h))                                      &
           + (sderk_bybarp-sderk_bybarm)*(dW*dW-h)/(4.0_wp*sqrt(h))  &
           + (sderk_bybarp-2.0_wp*sderk_b+sderk_bybarm)              &
             * (dW*h-I10)/(2.0_wp*h)                                 &
           + (sderk_bphip-sderk_bphim-sderk_bybarp+sderk_bybarm)     &
             * (dW*dW/3.0_wp-h)*dW/(4.0_wp*h)


  if ( (mf - 10*(mf/10)) .eq. 3 ) then
    sderk_c = sderk_c + 0.5_wp * sderk_a * h     
    call ord_2_iterator(t, h)
  else if ( (mf - 10*(mf/10)) .eq. 4 ) then
    call ord_4_iterator(t, h)
  else if ( (mf - 10*(mf/10)) .eq. 5 ) then
    call ord_6_iterator(t, h)
  end if

  return

end subroutine rk15_ito_imp_onestep


!----------------------------------------------------------------------
!
!  subroutine rk4_onestep(t, h)
!
!  Calculates one 4th-order deterministic Runge-Kutta iteration.
!  The parameter t is the time at the beginning of the interval of
!  length h. This subroutine requires the user-defined arrays
!  sderk_aphip, sderk_aphim, sderk_bphip, sderk_bphim, and sderk_ybar.
!
!-----------------------------------------------------------------------

subroutine rk4_onestep(t, h)

  implicit none

  real(wp), intent(in) :: t, h


  call sderk_func_a(t, sderk_y, sderk_aphip)
  sderk_ybar = sderk_y + 0.5_wp * h * sderk_aphip

  call sderk_func_a(t + 0.5_wp * h, sderk_ybar, sderk_aphim)
  sderk_ybar = sderk_y + 0.5_wp * h * sderk_aphim

  call sderk_func_a(t + 0.5_wp * h, sderk_ybar, sderk_bphip)
  sderk_ybar = sderk_y + h * sderk_bphip

  call sderk_func_a(t + h, sderk_ybar, sderk_bphim)
  sderk_y = sderk_y + (h/6.0_wp) * &
    (sderk_aphip + 2.0_wp * sderk_aphim + 2.0_wp * sderk_bphip + sderk_bphim)

  return

end subroutine rk4_onestep


!-----------------------------------------------------------------------
!
!  subroutine ord_2_iterator(t, h)
!
!  Performs iteration to convergence on the deterministic part of
!  the equations of motion for the implicit, order 2 deterministic
!  case.  Everything but the implicit part is stored in sderk_c,
!  and sderk_a stores the initial guess for the deterministic part.
!  Solution is by functional iteration; sderk_ybar contains the
!  updated solution and sderk_b is temporary storage to compute
!  residuals.
!
!-----------------------------------------------------------------------

subroutine ord_2_iterator(t, h)

  implicit none

  real(wp), intent(in) :: t, h

  integer  :: iter
  logical  :: stopflag
  real(wp) :: resid


  ! Make explicit guess for y_(n+1) and store as ybar
  sderk_ybar = sderk_y + 0.5_wp * sderk_a * h + sderk_c

  iter = 0
  stopflag = .false.
  do while ( .not. stopflag )
    iter = iter + 1
    call sderk_func_a(t + h, sderk_ybar, sderk_a)
    sderk_b = sderk_ybar
    sderk_ybar = sderk_y + 0.5_wp * sderk_a * h + sderk_c
    resid = sum(abs(sderk_b - sderk_ybar))

    if ( resid .lt. tol ) then
      stopflag = .true.
      if ( mf/10 .eq. 3 ) then
        write(0,*) 'Converged in ', iter, ' iterations at time t = ', t
      end if
    end if

    if ( ( .not. stopflag ) .and. ( iter .ge. maxit ) ) then
      write(0,*) 'Warning (SDERK): no convergence in ', iter, ' iterations'
      write(0,*) '  tol = ', tol, ', resid = ', resid
      stopflag = .true.
    end if      
  end do

  sderk_y = sderk_ybar

  return

end subroutine ord_2_iterator


!-----------------------------------------------------------------------
!
!  subroutine ord_4_iterator(t, h)
!
!  Performs iteration to convergence on the deterministic part of
!  the equations of motion for the implicit, order 4 deterministic      
!  case (see Iserles, A First Course in the Numerical Analysis of 
!  Differential Equations p. 47).  Everything but the implicit part 
!  is stored in sderk_c.  We will use sderk_a and sderk_b to store 
!  xi1 and xi2 arrays, and sderk_aphip, sderk_bphip, sderk_aphim, 
!  and sderk_bphim for temporary storage.
!  Solution is by functional iteration.
!
!-----------------------------------------------------------------------

subroutine ord_4_iterator(t, h)

  implicit none

  real(wp), intent(in) :: t, h

  integer  :: iter
  logical  :: stopflag
  real(wp) :: resid

  real(wp), dimension(2,2) :: ac
  real(wp), dimension(2)   :: bc, cc

  ac(:,1) = (/ 0.25_wp, 0.25_wp + sqrt(3.0_wp)/6.0_wp /)
  ac(:,2) = (/ 0.25_wp - sqrt(3.0_wp)/6.0_wp, 0.25_wp /)
  bc = 0.5_wp
  cc = (/ 0.5_wp - sqrt(3.0_wp)/6.0_wp, 0.5_wp + sqrt(3.0_wp)/6.0_wp /)


  ! initial guesses for xi1, xi2
  sderk_b = sderk_y + (sderk_a * h + sderk_c) * sum(ac(2,:))
  sderk_a = sderk_y + (sderk_a * h + sderk_c) * sum(ac(1,:))


  iter = 0
  stopflag = .false.
  do while ( .not. stopflag )
    iter = iter + 1
    call sderk_func_a(t + cc(1)*h, sderk_a, sderk_aphip)
    call sderk_func_a(t + cc(2)*h, sderk_b, sderk_bphip)
    sderk_aphim = sderk_a
    sderk_bphim = sderk_b
    sderk_a = sderk_y + h*(ac(1,1)*sderk_aphip+ac(1,2)*sderk_bphip) &
                        + sum(ac(1,:))*sderk_c
    sderk_b = sderk_y + h*(ac(2,1)*sderk_aphip+ac(2,2)*sderk_bphip) &
                        + sum(ac(2,:))*sderk_c
    resid = sum(abs(sderk_aphim-sderk_a)) + sum(abs(sderk_bphim-sderk_b))

    if ( resid .lt. tol ) then
      stopflag = .true.
      if ( mf/10 .eq. 3 ) then
        write(0,*) 'Converged in ', iter, ' iterations at time t = ', t
      end if
    end if

    if ( ( .not. stopflag ) .and. ( iter .ge. maxit ) ) then
      write(0,*) 'Warning (SDERK): no convergence in ', iter, ' iterations'
      write(0,*) '  tol = ', tol, ', resid = ', resid
      stopflag = .true.
    end if
  end do

  call sderk_func_a(t + cc(1)*h, sderk_a, sderk_aphip)
  call sderk_func_a(t + cc(2)*h, sderk_b, sderk_bphip)
  sderk_y = sderk_y + h*(bc(1)*sderk_aphip+bc(2)*sderk_bphip) + sderk_c

  return

end subroutine ord_4_iterator


!-----------------------------------------------------------------------
!
!  subroutine ord_6_iterator(t, h)
!
!  Performs iteration to convergence on the deterministic part of
!  the equations of motion for the implicit, order 6 deterministic
!  case (see Iserles, A First Course in the Numerical Analysis of
!  Differential Equations p. 47).  Everything but the implicit part
!  is stored in sderk_c.  We will use sderk_a, sderk_b, and sderk_ybar
!  to store the xi1, xi2, and xi3 arrays, and sderk_aybarp, 
!  sderk_bybarp, sderk_aphip, sderk_aybarm, sderk_bybarm, and
!  sderk_aphim for temporary storage.
!  Solution is by functional iteration.
!
!-----------------------------------------------------------------------

subroutine ord_6_iterator(t, h)

  implicit none

  real(wp), intent(in) :: t, h

  integer  :: iter
  logical  :: stopflag
  real(wp) :: resid

  real(wp), dimension(3,3) :: ac
  real(wp), dimension(3)   :: bc, cc

  ac(:,1) = (/ 5.0_wp/36.0_wp,                           & 
               5.0_wp/36.0_wp + sqrt(15.0_wp)/24.0_wp,   &
               5.0_wp/36.0_wp + sqrt(15.0_wp)/30.0_wp /)
  ac(:,2) = (/ 2.0_wp/9.0_wp - sqrt(15.0_wp)/15.0_wp,    &
               2.0_wp/9.0_wp,                            &
               2.0_wp/9.0_wp + sqrt(15.0_wp)/15.0_wp  /)
  ac(:,3) = (/ 5.0_wp/36.0_wp - sqrt(15.0_wp)/30.0_wp,   &
               5.0_wp/36.0_wp - sqrt(15.0_wp)/24.0_wp,   &
               5.0_wp/36.0_wp                         /)    
  bc = (/ 5.0_wp/18.0_wp, 4.0_wp/9.0_wp, 5.0_wp/18.0_wp /)     
  cc = (/ 0.5d0 - dsqrt(15.0d0)/10.0d0, 0.5d0,           &
          0.5d0 + dsqrt(15.0d0)/10.0d0                /)


  ! initial guesses for xi1, xi2, xi3
  sderk_ybar = sderk_y + (sderk_a * h + sderk_c) * sum(ac(3,:))
  sderk_b    = sderk_y + (sderk_a * h + sderk_c) * sum(ac(2,:))
  sderk_a    = sderk_y + (sderk_a * h + sderk_c) * sum(ac(1,:))


  iter = 0
  stopflag = .false.
  do while ( .not. stopflag )
    iter = iter + 1
    call sderk_func_a(t + cc(1)*h, sderk_a,    sderk_aybarp)
    call sderk_func_a(t + cc(2)*h, sderk_b,    sderk_bybarp)
    call sderk_func_a(t + cc(3)*h, sderk_ybar, sderk_aphip )
    sderk_aybarm = sderk_a
    sderk_bybarm = sderk_b
    sderk_aphim  = sderk_ybar
    sderk_a = sderk_y +                                              &
                h*(ac(1,1)*sderk_aybarp +                            &
                   ac(1,2)*sderk_bybarp +                            &
                   ac(1,3)*sderk_aphip   ) +                         &
                sum(ac(1,:)) * sderk_c
    sderk_b = sderk_y +                                              &
                h*(ac(2,1)*sderk_aybarp +                            &
                   ac(2,2)*sderk_bybarp +                            &
                   ac(2,3)*sderk_aphip   ) +                         &
                sum(ac(2,:)) * sderk_c
    sderk_ybar = sderk_y +                                           &
                h*(ac(3,1)*sderk_aybarp +                            &
                   ac(3,2)*sderk_bybarp +                            &
                   ac(3,3)*sderk_aphip   ) +                         &
                sum(ac(3,:)) * sderk_c
    resid = sum(abs(sderk_aybarm-sderk_a)) +                         &
            sum(abs(sderk_bybarm-sderk_b)) +                         &
            sum(abs(sderk_aphim-sderk_ybar))

    if ( resid .lt. tol ) then
      stopflag = .true.
      if ( mf/10 .eq. 3 ) then
        write(0,*) 'Converged in ', iter, ' iterations at time t = ', t
      end if
    end if

    if ( ( .not. stopflag ) .and. ( iter .ge. maxit ) ) then
      write(0,*) 'Warning (SDERK): no convergence in ', iter, ' iterations'
      write(0,*) '  tol = ', tol, ', resid = ', resid
      stopflag = .true.
    end if
  end do


  call sderk_func_a(t + cc(1)*h, sderk_a, sderk_aybarp)
  call sderk_func_a(t + cc(2)*h, sderk_b, sderk_bybarp)
  call sderk_func_a(t + cc(3)*h, sderk_ybar, sderk_aphip)
  sderk_y = sderk_y +                                                     &
            h*(bc(1)*sderk_aybarp+bc(2)*sderk_bybarp+bc(3)*sderk_aphip) + &
            sderk_c

  return

end subroutine ord_6_iterator


!-----------------------------------------------------------------------
!
!  function sderk_get_I1()
!
!  Returns a real vector containing the values of I1 used during the
!    last call to sderk.  The length is the value of nsteps used 
!    during the last call.
!
!-----------------------------------------------------------------------

function sderk_get_I1()
  implicit none
  real(wp), dimension(size(I1)) :: sderk_get_I1
  if ( .not. allocated(I1) ) then
    write(0,*) 'Error (SDERK_GET_I1): SDERK not yet called'
    stop
  end if
  sderk_get_I1 = I1
end function sderk_get_I1


!-----------------------------------------------------------------------
!
!  function sderk_get_I10()     
!
!  Returns a real vector containing the values of I10 used during the
!    last call to sderk.  The length is the value of nsteps used 
!    during the last call.
!
!-----------------------------------------------------------------------

function sderk_get_I10()
  implicit none
  real(wp), dimension(size(I10)) :: sderk_get_I10
  if ( .not. allocated(I10) ) then
    write(0,*) 'Error (SDERK_GET_I10): SDERK not yet called'
    stop
  end if
  sderk_get_I10 = I10
end function sderk_get_I10

  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  Now comes an embedded version of the random_pl random number
!  generator (multi routines only)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!=======================================================================
!
!  subroutine init_rand_pl_multi(m, optional seed1, optional seed2,
!                                optional seed3)
!
!  I: m          integer      selects which generator to initialize
!     seed1,2,3  int(rpk)     optional seeds for initializing generator
!
!  O: globals only
!
!  Analogous routine to init_rand_pl, to initialize one of the set
!  of multiple random number generators, specified by m.
!
!=======================================================================

subroutine init_rand_pl_multi(m, seed1, seed2, seed3)

  implicit none

  integer, intent(in) :: m
  integer(rpk), optional, intent(in) :: seed1, seed2, seed3
  integer(rpk) :: lseed1, lseed2, lseed3, firstseed

  ! set seeds
  if ( present(seed1) .or. present(seed2) .or. present(seed3) ) then
    if ( present(seed3) ) firstseed = seed3
    if ( present(seed2) ) firstseed = seed2
    if ( present(seed1) ) firstseed = seed1
    lseed1 = firstseed; lseed2 = firstseed; lseed3 = firstseed
    if ( present(seed1) ) lseed1 = seed1
    if ( present(seed2) ) lseed2 = seed2
    if ( present(seed3) ) lseed3 = seed3
  else
    lseed1 = dseed1+m; lseed2 = dseed2+m; lseed3 = dseed3+m
  end if

  call associate_multi(m)
  call init_rand_pl_main(lseed1, lseed2, lseed3)

end subroutine init_rand_pl_multi


!=======================================================================
!
!  subroutine allocate_rand_pl_multi(m, optional seed1, optional seed2,
!                              optional seed3)
!
!  I: m          integer      selects how many generators to allocate
!     seed1,2,3  int(rpk)     optional seeds for initializing generator
!
!  O: globals only
!
!  Allocate storage for m multiple random number generators.  Also
!  initializes all generators, as follows: 
!  the multi generators 1 through m are set using the specified seeds
!  incremented by the generator number.
!
!  Any existing set of multiple generators will be obliterated.
!  
!=======================================================================

subroutine allocate_rand_pl_multi(m, seed1, seed2, seed3)

  implicit none

  integer, intent(in) :: m
  integer(rpk), optional, intent(in) :: seed1, seed2, seed3
  integer(rpk) :: lseed1, lseed2, lseed3, firstseed
  integer :: j
  integer(rpk) :: jrpk
  
  ! set seeds
  if ( present(seed1) .or. present(seed2) .or. present(seed3) ) then
    if ( present(seed3) ) firstseed = seed3
    if ( present(seed2) ) firstseed = seed2
    if ( present(seed1) ) firstseed = seed1
    lseed1 = firstseed; lseed2 = firstseed; lseed3 = firstseed
    if ( present(seed1) ) lseed1 = seed1
    if ( present(seed2) ) lseed2 = seed2
    if ( present(seed3) ) lseed3 = seed3
  else 
    lseed1 = dseed1; lseed2 = dseed2; lseed3 = dseed3
  end if 
  
  if ( allocated( initialized_m ) ) then
    deallocate( lecuyer_stor_m, scramble_stor_m, fib_stor_m, mersenne_stor_m)
    deallocate( initialized_m, fib_ctr_m, bd_lastout_m, mersenne_ptr_m )
  end if
  allocate( lecuyer_stor_m(lecuyer_stor_sz, m) )
  allocate( scramble_stor_m(scramble_stor_sz, m) )
  allocate( fib_stor_m(fib_stor_sz, m), mersenne_stor_m(mersenne_stor_sz, m) )
  allocate( initialized_m(m), fib_ctr_m(m) )
  allocate( bd_lastout_m(m), mersenne_ptr_m(m) )

  do j = 1, m
    jrpk = j
    call associate_multi(j)
    call init_rand_pl_main(lseed1+jrpk, lseed2+jrpk, lseed3+jrpk)
    initialized = 1
  end do

end subroutine allocate_rand_pl_multi


!=======================================================================
!
!  real(kind=wp) function rand_pl_multi(m)
!
!  I: m              integer          selects which generator to use
!  
!  O: rand_pl_multi  real(kind=wp)    random number in (0,1)
!
!  This is a wrapper for rand_pl_main that associates pointers for
!  the multiple generators before calling the main routine. 
!
!=======================================================================

function rand_pl_multi(m) 
  implicit none
  integer, intent(in) :: m
  real(wp) :: rand_pl_multi
  if ( curr_ptr .ne. m ) call associate_multi(m)
  if ( initialized .ne. 1 ) call init_rand_pl_multi(m)
  rand_pl_multi = rand_pl_main()
end function rand_pl_multi


!=======================================================================
!
!  subroutine associate_multi(m)
!
!  I: m        integer        selects which of the set of multiple 
!                             generators to use
!
!  O: globals only
!
!  Associates pointers so that operations are on one of the multiple
!  generators, specified by the integer n.
!  
!=======================================================================
  
subroutine associate_multi(m)

  implicit none

  integer, intent(in) :: m

  if ( .not. allocated( initialized_m ) ) then
    write(0,*) 'Error (ASSOCIATE_MULTI): attempt to use multiple ', &
               'generators before initialization'
    stop
  end if

  if ( m .lt. 1 .or. m .gt. size(initialized_m) ) then
    write(0,*) 'Error (ASSOCIATE_MULTI): selector is out of allocated range'
    write(0,*) '  selector = ', m, ', max = ', size(initialized_m)
    stop
  end if
  
  curr_ptr = m
  lecuyer_stor  => lecuyer_stor_m(:,m)
  scramble_stor => scramble_stor_m(:,m)
  fib_stor      => fib_stor_m(:,m)
  mersenne_stor => mersenne_stor_m(:,m)
  initialized   => initialized_m(m)
  fib_ctr       => fib_ctr_m(m)
  bd_lastout    => bd_lastout_m(m)
  mersenne_ptr  => mersenne_ptr_m(m)

end subroutine associate_multi


!=======================================================================
!
!  integer(kind=rpk) function mersenne()
!
!  I: only through module globals
!
!  O: mersenne  int(rpk)    random integer, in the range
!                           [0, 2**31 - 1] inclusive
!
!  Implements the Mersenne twister generator, returning one random
!  integer.
!
!=======================================================================

function mersenne()

  implicit none

  integer(rpk) :: mersenne

  integer(rpk), parameter :: n = 624, m = 397
  integer(rpk), parameter :: unity      =           1_rpk  ! 00000001
  integer(rpk), parameter :: matrix_a   = -1727483681_rpk  ! 9908b0df
  integer(rpk), parameter :: upp_p1     = -2147483647_rpk  ! 80000001
  integer(rpk), parameter :: upp_mask   =    upp_p1-unity  ! 80000000
  integer(rpk), parameter :: low_mask   =  2147483647_rpk  ! 7fffffff
  integer(rpk), parameter :: tmp_mask_b = -1658038656_rpk  ! 9d2c5680
  integer(rpk), parameter :: tmp_mask_c =  -272236544_rpk  ! efc60000
  integer(rpk), dimension(0:1), parameter :: mag01 = (/ 0_rpk, matrix_a /)

  integer, parameter :: shift_u = -11
  integer, parameter :: shift_s =   7
  integer, parameter :: shift_t =  15
  integer, parameter :: shift_l = -18

  integer(rpk) :: y
  integer :: k

  if ( mersenne_ptr .ge. n+1 ) then
    do k = 1, n-m
      y = ior( iand(mersenne_stor(k),   upp_mask), &
               iand(mersenne_stor(k+1), low_mask) )
      mersenne_stor(k) = ieor( ieor( mersenne_stor(k+m), ishft(y,-1) ), &
                               mag01( iand(y, unity) ) )
    end do

    do k = n-m+1, n-1
      y = ior( iand(mersenne_stor(k),   upp_mask), &
               iand(mersenne_stor(k+1), low_mask) )
      mersenne_stor(k) = ieor( ieor( mersenne_stor(k+m-n), ishft(y,-1) ), &
                               mag01( iand(y, unity) ) )
    end do

    y = ior( iand(mersenne_stor(n), upp_mask), &
             iand(mersenne_stor(1), low_mask) )
    mersenne_stor(n) = ieor( ieor( mersenne_stor(m), ishft(y,-1) ), &
                             mag01( iand(y, unity) ) )

    mersenne_ptr = 1
  end if

  y = mersenne_stor(mersenne_ptr)
  mersenne_ptr = mersenne_ptr + 1
  y = ieor( y, ishft(y, shift_u) )
  y = ieor( y, iand(ishft(y, shift_s), tmp_mask_b) )
  y = ieor( y, iand(ishft(y, shift_t), tmp_mask_c) )
  y = ieor( y, ishft(y, shift_l) )

  ! shift off 1 bit, so we get a 31-bit integer
  mersenne = ishft(y, -1)

end function mersenne


!=======================================================================
!
!  subroutine init_mersenne(seed)
!
!  I: seed     int(rpk)       initialization seed for generator,
!                              any nonzero integer is ok
!
!  O: none
!
!  Initializes Mersenne twister generator as in the sample code by
!  Matsumoto and Nishimura, using Knuth's LCG.
!
!=======================================================================
  
subroutine init_mersenne(seed)
  
  implicit none 
  
  integer(rpk), intent(in) :: seed

  integer(rpk), parameter :: mask32 = -1_rpk ! ffffffff

  ! check seed
  if ( iand(seed, mask32) .eq. 0_rpk ) then
    write(0,*) 'Error (INIT_MERSENNE): nonzero seed required'
    stop
  end if

  ! just in case the storage uses more than 32 bits
  mersenne_stor(1) = iand( seed, mask32 )
  
  do mersenne_ptr = 2, mersenne_stor_sz
    mersenne_stor(mersenne_ptr) = 1812433253_rpk * &
      ieor( mersenne_stor(mersenne_ptr-1), &
            ishft( mersenne_stor(mersenne_ptr-1), -30 ) ) + mersenne_ptr - 1
    mersenne_stor(mersenne_ptr) = iand( mersenne_stor(mersenne_ptr), mask32 )
  end do

  mersenne_ptr = mersenne_stor_sz + 1

end subroutine init_mersenne
  

!=======================================================================
!  
!  integer(kind=rpk) function lecuyer()
!
!  I: only through module globals
!
!  O: lecuyer   int(rpk)    random integer, in the range
!                           [0, 2**31 - 2] inclusive
!
!  Implements one iteration of the L'Ecuyer generator.
!
!=======================================================================

function lecuyer()

  implicit none

  integer(rpk) :: lecuyer

  integer(rpk), parameter :: m1  = 2147483647_rpk ! 2**31 - 1
  integer(rpk), parameter :: m2  = 2145483479_rpk

  integer(rpk), parameter :: a12 =      63308_rpk
  integer(rpk), parameter :: a13 =    -183326_rpk
  integer(rpk), parameter :: a21 =      86098_rpk
  integer(rpk), parameter :: a23 =    -539608_rpk

  integer(rpk), parameter :: q12 =      33921_rpk
  integer(rpk), parameter :: q13 =      11714_rpk
  integer(rpk), parameter :: q21 =      24919_rpk
  integer(rpk), parameter :: q23 =       3976_rpk

  integer(rpk), parameter :: r12 =      12979_rpk
  integer(rpk), parameter :: r13 =       2883_rpk
  integer(rpk), parameter :: r21 =       7417_rpk
  integer(rpk), parameter :: r23 =       2071_rpk

  integer(rpk) :: h, p12, p13, p21, p23
  integer(rpk), pointer :: x10, x11, x12, x20, x21, x22
  x10 => lecuyer_stor(1); x11 => lecuyer_stor(2); x12 => lecuyer_stor(3)
  x20 => lecuyer_stor(4); x21 => lecuyer_stor(5); x22 => lecuyer_stor(6)

  ! Component 1
  h = x10 / q13
  p13 = -a13 * (x10 - h * q13) - h * r13
  h = x11 / q12
  p12 =  a12 * (x11 - h * q12) - h * r12
  if ( p13 .lt. 0 ) p13 = p13 + m1
  if ( p12 .lt. 0 ) p12 = p12 + m1
  x10 = x11
  x11 = x12
  x12 = p12 - p13
  if ( x12 .lt. 0 ) x12 = x12 + m1

  ! Component 2
  h = x20 / q23
  p23 = -a23 * (x20 - h * q23) - h * r23
  h = x22 / q21
  p21 =  a21 * (x22 - h * q21) - h * r21
  if ( p23 .lt. 0 ) p23 = p23 + m2
  if ( p21 .lt. 0 ) p21 = p21 + m2
  x20 = x21
  x21 = x22
  x22 = p21 - p23
  if ( x22 .lt. 0 ) x22 = x22 + m2

  ! Combination
  if ( x12 .le. x22 ) then
    lecuyer = x12 - x22 + m1 - 1
  else
    lecuyer = x12 - x22 - 1
  end if

end function lecuyer


!=======================================================================
!
!  subroutine init_lecuyer(seed)
!
!  I: seed     int(rpk)       initialization seed for generator,
!                              different seeds give distinct sequences
!                              by partitioning the full sequence;
!                              any number in [0, 2**31 - 2] is ok
!
!  O: none
!
!  Initializes lecuyer generator by setting the seeds from a single
!  integer, by partitioning the full sequence and jumping ahead
!  by an amount proportional to the seed.
!
!=======================================================================

subroutine init_lecuyer(seed)

  implicit none
  
  integer(rpk), intent(in) :: seed

  integer(rpk) :: lseed

  integer(rpk), parameter :: m1  = 2147483647_rpk ! 2**31 - 1
  integer(rpk), parameter :: m2  = 2145483479_rpk

  integer(rpk), parameter :: a12 =      63308_rpk
  integer(rpk), parameter :: a13 =    -183326_rpk + m1
  integer(rpk), parameter :: a21 =      86098_rpk
  integer(rpk), parameter :: a23 =    -539608_rpk + m2

  integer(rpk), parameter :: logk = 184 - 32
  integer(rpk), dimension(3,3) :: A1, A2, A1k, A2k, A1ks, A2ks
  integer(rpk), dimension(3) :: tmparr

  integer :: j

  
  ! check seed
  if ( seed .lt. 0_rpk .or. seed .gt. 2147483646_rpk ) then
    write(0,*) 'Error (INIT_LECUYER): seed outside acceptable range'
    stop
  end if
  lseed = seed + 1

  ! single-iteration matrices
  A1(1,:) = (/   0_rpk,   1_rpk,   0_rpk /)
  A1(2,:) = (/   0_rpk,   0_rpk,   1_rpk /)
  A1(3,:) = (/   0_rpk,     a12,     a13 /)
  A2(1,:) = (/   0_rpk,   1_rpk,   0_rpk /)
  A2(2,:) = (/   0_rpk,   0_rpk,   1_rpk /)
  A2(3,:) = (/     a21,   0_rpk,     a23 /)

  ! compute k-iteration matrices, to jump between different partitions
  A1k = A1
  A2k = A2
  do j = 1, logk
    A1k = modmatsq(A1k, m1)  ! A1k = A1k * A1k (mod m1)
    A2k = modmatsq(A2k, m2)  ! A2k = A2k * A2k (mod m2)
  end do

  ! compute (A**k)**lseed matrices, to jump ahead to the location 
  !   marked by lseed
  A1ks = divconq(A1k, lseed, m1)
  A2ks = divconq(A2k, lseed, m2)

  ! set default seeds
  lecuyer_stor = (/ 1852689663, 1962642687, 580869375, &
                    2039711750, 1671394257, 879888250  /)

  ! advance by k*seed iterations (tmparr used to work around ifc bug)
  tmparr = lecuyer_stor(1:3)
  lecuyer_stor(1:3) = modmvmul(A1ks, tmparr, m1)
  tmparr = lecuyer_stor(4:6)
  lecuyer_stor(4:6) = modmvmul(A2ks, tmparr, m2)

end subroutine init_lecuyer


  !!!!
  !!!! routines for integer modular matrix/vector/scalar arithmetic
  !!!!

  ! divide-and-conquer to calculate A**j mod m
  recursive function divconq(A, j, m) result (Aj)
    implicit none
    integer(rpk), dimension(:,:), intent(in) :: A
    integer(rpk), intent(in) :: j, m
    integer(rpk), dimension(size(A,1),size(A,2)) :: Aj
    if ( j .eq. 1 ) then
      Aj = A
    else if ( isodd(j) ) then
      Aj = modmatmul(A, divconq(A, j-1_rpk, m), m)
    else
      Aj = modmatsq(divconq(A, j/2_rpk, m), m)
    end if
  end function divconq
    
  ! check if integer scalar is odd
  function isodd(a)
    implicit none
    integer(rpk), intent(in) :: a
    logical :: isodd
    isodd = (iand(a, 1_rpk) .eq. 1_rpk)
  end function isodd

  ! wrapper for modmatmul to square a matrix
  function modmatsq(A, m)
    implicit none
    integer(rpk), dimension(:,:), intent(in) :: A
    integer(rpk), intent(in) :: m
    integer(rpk), dimension(size(A,1),size(A,2)) :: modmatsq
    modmatsq = modmatmul(A, A, m)
  end function modmatsq

  ! calculates the matrix product A * B mod m, assuming nonnegative
  ! entries; A and B are assumed square
  function modmatmul(A, B, m)
    implicit none
    integer(rpk), dimension(:,:), intent(in) :: A, B
    integer(rpk), intent(in) :: m
    integer(rpk), dimension(size(A,1),size(A,2)) :: modmatmul
    integer(rpk) :: n, j
    n = size(A,1)
    if ( n.ne.size(A,2) .or. n.ne.size(B,1) .or. n.ne.size(B,2) ) then
      write(0,*) 'Error (MODMATMUL): nonsquare or incompatible matrix sizes'
      stop
    end if
    do j = 1, n
      modmatmul(:,j) = modmvmul(A, B(:,j), m)
    end do
  end function modmatmul

  ! calculates the matrix-vector product A * x mod m, assuming 
  ! nonnegative entries; A is assumed square
  function modmvmul(A, x, m)
    implicit none
    integer(rpk), dimension(:,:), intent(in) :: A
    integer(rpk), dimension(:),   intent(in) :: x
    integer(rpk), intent(in) :: m
    integer(rpk), dimension(size(A,1)) :: modmvmul
    integer(rpk) :: n, j
    n = size(A,1)
    if ( n.ne.size(A,2) .or. n.ne.size(x) ) then
      write(0,*) 'Error (MODMVMUL): nonsquare matrix or incompatible vector'
      stop
    end if
    do j = 1, n
      modmvmul(j) = moddotprod(A(j,:), x, m)
    end do
  end function modmvmul

  ! calculates the vector dot product x * y mod m, assuming 
  ! nonnegative entries
  function moddotprod(x, y, m)
    implicit none
    integer(rpk), dimension(:),   intent(in) :: x, y
    integer(rpk), intent(in) :: m
    integer(rpk) :: moddotprod
    integer(rpk) :: n, j
    n = size(x)
    if ( n.ne.size(y) ) then
      write(0,*) 'Error (MODDOTPROD): incompatible vector sizes'
      stop
    end if
    moddotprod = modmult(x(1), y(1), m)
    do j = 2, n
      moddotprod = modadd( moddotprod, modmult(x(j), y(j), m), m )
    end do
  end function moddotprod

  ! calculates x+y mod m without overflow, assuming x, y >= 0
  function modadd(x, y, m)
    implicit none
    integer(rpk), intent(in) :: x, y, m
    integer(rpk) :: modadd
    modadd = (x - m) + y; if ( modadd .lt. 0 ) modadd = modadd + m
  end function modadd
    
  ! calculates x*y mod m without overflow via the decomposition
  ! algorithm; see L'Ecuyer and Cote, ACM Trans. Math. Soft. 17, 98 
  ! (1991). Assumes x, y >= 0.  Note the redundant p = p+m statements
  ! that were necessary to work around an intel compiler bug.
  function modmult(x, y, m)
    implicit none
    integer(rpk), intent(in) :: x, y, m
    integer(rpk) :: modmult
    integer(rpk), parameter :: h = 32768 ! for 32-bit integers
    integer(rpk) :: a, s, a0, a1, q, qh, rh, k, p
    if ( x .lt. 0 .or. y .lt. 0 ) then
      write(0,*) 'Error (MODMULT): inputs must be positive'
      stop
    end if
    if ( x .eq. 0 .or. y .eq. 0 ) then
      modmult = 0; return
    end if
    a = x; s = y
    if ( a .lt. h ) then
      a0 = a; p = 0
    else
      a1 = a / h; a0 = a - h*a1
      qh = m / h; rh = m - h*qh
      if ( a1 .ge. h ) then
        a1 = a1 - h; k = s / qh
        p = h * (s - k*qh) - k * rh
        if ( p .lt. 0 ) p = p + m
        do while ( p .lt. 0 ); p = p + m; end do
      else
        p = 0
      end if
      if ( a1 .ne. 0 ) then
        q = m / a1; k = s / q
        p = p - k * (m - a1*q)
        if ( p .gt. 0 ) p = p - m
        p = p + a1 * (s - k*q)
        if ( p .lt. 0 ) p = p + m
        do while ( p .lt. 0 ); p = p + m; end do
      end if
      k = p / qh; p = h * (p - k*qh) - k * rh
      if ( p .lt. 0 ) p = p + m
      do while ( p .lt. 0 ); p = p + m; end do
    end if
    if ( a0 .ne. 0 ) then
      q = m / a0; k = s / q
      p = p - k * (m - a0*q)
      if ( p .gt. 0 ) p = p - m
      p = p + a0 * (s - k*q)
      if ( p .lt. 0 ) p = p + m
      do while ( p .lt. 0 ); p = p + m; end do
    end if
    modmult = p
  end function modmult


!=======================================================================
!  
!  integer(kind=rpk) function fibonacci(optional force_init)
!
!  I: force_init  logical    (optional) force a new set of 1009 
!                            numbers to be generated
!
!  O: fibonacci   int(rpk)    random integer, in the range
!                           [0, 2**31 - 2] (note that the lsb is 
!                                           always 0)
!
!  Implements one iteration of Knuth's subtractive lagged Fibonacci
!    generator, with skipping to improve randomness.
!
!=======================================================================

function fibonacci(force_init)

  implicit none

  logical, optional :: force_init

  integer(rpk) :: fibonacci

  integer(rpk), parameter :: lag = 37, sbst = 100, total = 1009
  integer(rpk), parameter :: mdls = 1073741824_rpk ! 2**30
  logical :: reinit
  integer :: j

  ! check for forcing reinitialization
  reinit = .false.
  if ( present( force_init ) ) reinit = force_init

  ! generate 1009 numbers, if needed
  if ( fib_ctr .gt. sbst .or. reinit ) then
    aa(1:sbst) = fib_stor(1:sbst)
    do j = sbst+1, total
      aa(j) = iand(aa(j-sbst) - aa(j-lag), mdls-1_rpk)
    end do
    do j = 1, lag
      fib_stor(j) = iand(aa(total+j-sbst) - aa(total+j-lag), mdls-1_rpk)
    end do
    do j = lag+1, sbst
      fib_stor(j) = iand(aa(total+j-sbst) - fib_stor(j-lag), mdls-1_rpk)
    end do
    fib_ctr = 1
  end if

  fibonacci = 2 * fib_stor(fib_ctr)
  fib_ctr = fib_ctr + 1

end function fibonacci
  

!=======================================================================
!  
!  subroutine init_fibonacci(seed)
!
!  I: seed    int(rpk)   initialization seed, in the range
!                        [0, 2**30 - 3 = 1 073 741 821]
!
!  O: only through module globals
!
!  Initializes Knuth's Fibonacci generator, according to his
!    recommendations.
!
!=======================================================================

subroutine init_fibonacci(seed)

  implicit none

  integer(rpk), intent(in) :: seed

  integer(rpk), parameter :: tt = 70, lag = 37 
  integer(rpk), parameter :: mdls = 1073741824_rpk ! 2**30
  integer(rpk), parameter :: maxseed = mdls-3_rpk
  integer, parameter :: kk = fib_stor_sz, kkk = kk+kk-1
  integer(rpk), dimension(kkk) :: xtmp
  integer :: t, j, ss
  integer(rpk) :: dummy


  ! check seed
  if ( seed .lt. 1 .or. seed .gt. maxseed ) then
    write(0,*) 'Error (INIT_FIBONACCI): inappropriate seed'
  end if

  ! set pointer
  fib_ctr = 1

  ! do knuthian stuff
  ss = seed - mod(seed,2_rpk) + 2
  do j = 1, kk
    xtmp(j) = ss
    ss = ss + ss
    if ( ss .ge. mdls ) ss = ss - mdls + 2
  end do
  xtmp(2) = xtmp(2) + 1
  ss = seed
  t = tt - 1
10 continue
  do j = kk, 2, -1
    xtmp(j+j-1) = xtmp(j)
    xtmp(j+j-2) = 0
  end do
  do j = kkk, kk+1, -1
    xtmp(j-(kk-lag)) = xtmp(j-(kk-lag)) - xtmp(j)
    if ( xtmp(j-(kk-lag) ) .lt. 0)  xtmp(j-(kk-lag)) = xtmp(j-(kk-lag)) + mdls
    xtmp(j-kk) = xtmp(j-kk) - xtmp(j)
    if ( xtmp(j-kk) .lt. 0 )  xtmp(j-kk) = xtmp(j-kk) + mdls
  end do
  if ( mod(ss, 2) .eq. 1 ) then
    do j = kk, 1, -1
      xtmp(j+1) = xtmp(j)
    end do
    xtmp(1) = xtmp(kk+1)
    xtmp(lag+1) = xtmp(lag+1) - xtmp(kk+1)
    if ( xtmp(lag+1) .lt. 0 )  xtmp(lag+1) = xtmp(lag+1) + mdls
  end if
  if ( ss .ne. 0 ) then
    ss = ss / 2
  else
    t = t - 1
  end if
  if ( t .gt. 0 ) go to 10
  do j = 1, lag
    fib_stor(j+kk-lag) = xtmp(j)
  end do
  do j = lag+1, kk
    fib_stor(j-lag) = xtmp(j)
  end do

  ! discard the first 2020 or so numbers
  do j = 1, 2
    dummy = fibonacci(force_init=.true.)
  end do

end subroutine init_fibonacci


!=======================================================================
!
!  subroutine init_rand_pl_main(seed1, seed2, seed3)
!
!  I: seed1,2,3  int(rpk)     seeds for initializing generator
!
!  O: through module globals
!
!  Initializes all 3 random number generators, with the
!  specified seeds (1 for Fibonacci, 2 for L'Ecuyer, 3 for Mersenne).  
!
!=======================================================================

subroutine init_rand_pl_main(seed1, seed2, seed3)

  implicit none

  integer(rpk), intent(in) :: seed1, seed2, seed3
  integer :: j


  ! check method flag
  if ( rand_pl_mf_pri .lt. 1 .or. rand_pl_mf_pri .gt. 3 .or. &
       rand_pl_mf_mix .lt. 0 .or. rand_pl_mf_mix .gt. 3 .or. &
       rand_pl_mf_scr .lt. 0 .or. rand_pl_mf_scr .gt. 4 .or. &
       rand_pl_mf_pri .eq. rand_pl_mf_mix .or. &
       rand_pl_mf_pri .eq. rand_pl_mf_scr .or. &
       (rand_pl_mf_mix .ne. 0 .and. rand_pl_mf_mix .eq. rand_pl_mf_scr) ) then
    write(0,*) 'Error (INIT_RAND_PL): illegal rand_pl_mf'
    stop
  end if

  ! initialize all generators
  call init_fibonacci(seed1)
  call init_lecuyer(seed2)
  call init_mersenne(seed3)

  ! fill scramble table
  bd_lastout = 0
  scramble_stor = 0
  if ( rand_pl_mf_scr .ne. 0 ) then
    do j = 1, scramble_stor_sz
      scramble_stor(j) = rand_pl_raw()
    end do
    if ( rand_pl_mf_scr .eq. 4 ) bd_lastout = rand_pl_raw()
  end if

  initialized = 1

end subroutine init_rand_pl_main


!=======================================================================
!
!  real(kind=wp) function rand_pl_raw()
!
!  I: none
!
!  O: rand_pl     integer(rpk)         random integer 
!                              
!  Constructs a random integer from the proper generator(s) and
!  subtractively combines them, if appropriate.  Basically this
!  does everything short of shuffling the output numbers.
!
!=======================================================================

function rand_pl_raw() result(r)

  implicit none

  integer(rpk) :: r

  integer(rpk), parameter :: modmask = 2147483647_rpk  ! 7fffffff = 2**31-1


  if ( rand_pl_mf_pri .eq. 1 ) then
    r = fibonacci()
  else if ( rand_pl_mf_pri .eq. 2 ) then
    r = lecuyer()
  else if ( rand_pl_mf_pri .eq. 3 ) then
    r = mersenne()
  end if

  if ( rand_pl_mf_mix .ne. 0 ) then
    if ( rand_pl_mf_mix .eq. 1 ) then
      r = iand( r - fibonacci(), modmask )
    else if ( rand_pl_mf_mix .eq. 2 ) then
      r = iand( r - lecuyer(), modmask )
    else if ( rand_pl_mf_mix .eq. 3 ) then
      r = iand( r - mersenne(), modmask )
    end if
  end if

end function rand_pl_raw


!=======================================================================
!
!  real(kind=wp) function rand_pl_main()
!
!  I: none
!
!  O: rand_pl     real(kind=wp)         random number in (0,1)
!
!=======================================================================

function rand_pl_main()

  implicit none

  real(wp) :: rand_pl_main

  integer(rpk), parameter :: fibonacci_range = 1073741823_rpk ! 2**31 - 2
  real(wp), parameter :: fibonacci_mult = 4.65661287307739e-10_wp ! 1/(2**31)
  integer(rpk), parameter :: lecuyer_range   = 2147483646_rpk ! 2**31 - 2
  real(wp), parameter :: lecuyer_mult = 4.65661287307739e-10_wp ! 1/(2**31)
  integer(rpk), parameter :: mersenne_range  = 2147483647_rpk ! 2**31 - 1
  real(wp), parameter :: mersenne_mult = 4.65661287090899e-10_wp ! 1/(2**31 + 1)

  real(wp), parameter :: gen_mult = 4.65661287090899e-10_wp ! 1/(2**31 + 1)

  integer :: indx

 
  if ( rand_pl_mf_scr .eq. 0 ) then
    rand_pl_main = (real(rand_pl_raw(), kind=wp) + 1) * gen_mult

  else if ( rand_pl_mf_scr .eq. 1 ) then
    ! scramble output using Fibonacci output
    indx = ceiling(fibonacci() * gen_mult * scramble_stor_sz)
    rand_pl_main = (real(scramble_stor(indx), kind=wp) + 1) * gen_mult
    scramble_stor(indx) = rand_pl_raw()

  else if ( rand_pl_mf_scr .eq. 2 ) then
    ! scramble output using L'Ecuyer output
    indx = ceiling(lecuyer() * gen_mult * scramble_stor_sz)
    rand_pl_main = (real(scramble_stor(indx), kind=wp) + 1) * gen_mult
    scramble_stor(indx) = rand_pl_raw()

  else if ( rand_pl_mf_scr .eq. 3 ) then
    ! scramble output using Mersenne output
    indx = ceiling(mersenne() * gen_mult * scramble_stor_sz)
    rand_pl_main = (real(scramble_stor(indx), kind=wp) + 1) * gen_mult
    scramble_stor(indx) = rand_pl_raw()

  else if ( rand_pl_mf_scr .eq. 4 ) then
    ! scramble output using Bays-Durham scrambling
    indx = ceiling(bd_lastout * gen_mult * scramble_stor_sz)
    bd_lastout = scramble_stor(indx)
    rand_pl_main = (real(bd_lastout, kind=wp) + 1) * gen_mult
    scramble_stor(indx) = rand_pl_raw()

  end if

end function rand_pl_main


!=======================================================================
!
!  real(kind=wp) function rand_pl_vec_multi(m, n)
!
!  I: m            default integer      selects which generator to use
!     n            default integer      length of output vector
!
!  O: rand_pl_vec  real(kind=wp) array  array of random numbers in (0,1)
!
!  Single wrapper for rand_pl_vec_main()
!
!=======================================================================

function rand_pl_vec_multi(m, n)
  implicit none
  integer, intent(in) :: m, n
  real(wp), dimension(n) :: rand_pl_vec_multi
  if ( curr_ptr .ne. m ) call associate_multi(m)
  rand_pl_vec_multi = rand_pl_vec_main(n)
end function rand_pl_vec_multi


!=======================================================================
!
!  real(kind=wp) function rand_pl_vec_main(n)
!
!  I: n            default integer      length of output vector
!
!  O: rand_pl_vec_main  real(kind=wp) array  array of random numbers 
!                                            in (0,1)
!
!=======================================================================

function rand_pl_vec_main(n)
  implicit none
  integer, intent(in) :: n
  real(wp), dimension(n) :: rand_pl_vec_main
  integer :: j
  do j = 1, n
    rand_pl_vec_main(j) = rand_pl_main()
  end do
end function rand_pl_vec_main


!=======================================================================
!
!  real(kind=wp) function nrand_pl_vec_multi(m, n)
!
!  I: m            default integer      selects which generator to use
!     n            default integer      length of output vector, must
!                                         be even
!
!  O: nrand_pl_vec real(kind=wp) array  array of standard normal
!                                         deviates
!  
!  Multiple wrapper for nrand_pl_vec_main
!  
!=======================================================================

function nrand_pl_vec_multi(m, n)
  implicit none
  integer, intent(in) :: m, n
  real(wp), dimension(n) :: nrand_pl_vec_multi
  if ( curr_ptr .ne. m ) call associate_multi(m)
  nrand_pl_vec_multi = nrand_pl_vec_main(n)
end function nrand_pl_vec_multi


!=======================================================================
!
!  real(kind=wp) function nrand_pl_vec_main(n)
!
!  I: n            default integer      length of output vector, must
!                                         be even
!
!  O: nrand_pl_vec_main real(kind=wp) array  array of standard normal
!                                            deviates
!
!=======================================================================

function nrand_pl_vec_main(n) result(nd)

  implicit none

  integer, intent(in) :: n
  real(wp), dimension(n) :: nd

  real(wp), parameter :: pi = 3.141592653589793_wp, twopi = 2*pi
  integer :: j
  real(wp) :: tmp

  
  if ( (n/2)*2 .ne. n ) then
    write(0,*) 'Error (NRAND_PL_VEC): requested length not even'
    stop
  end if

  nd = rand_pl_vec_main(n)

  do j = 1, n, 2
    tmp = sqrt((-2.0_wp) * log(nd(j)))
    nd(j)   = cos(twopi * nd(j+1)) * tmp
    nd(j+1) = sin(twopi * nd(j+1)) * tmp
  end do

end function nrand_pl_vec_main


!=======================================================================
!
!  subroutine get_rand_pl_state_multi(m, ivec)
!
!  I: m      integer            selects which generator to use
!
!  O: ivec   int(rpk) array     contains sufficient information to
!                               restore the state of the generator
!
!  Single wrapper for get_rand_pl_state_main
!
!=======================================================================
    
subroutine get_rand_pl_state_multi(m, ivec)
  implicit none
  integer(rpk), dimension(rand_pl_state_sz), intent(out) :: ivec
  integer, intent(in) :: m
  if ( curr_ptr .ne. m ) call associate_multi(m)
  call get_rand_pl_state_main(ivec)
end subroutine get_rand_pl_state_multi


!=======================================================================
!
!  subroutine get_rand_pl_state_main(ivec)
!
!  I: only through globals
!
!  O: ivec   int(rpk) array     contains sufficient information to
!                               restore the state of the generator
!
!=======================================================================

subroutine get_rand_pl_state_main(ivec)

  implicit none

  integer(rpk), dimension(rand_pl_state_sz), intent(out) :: ivec
  integer :: ptr1, ptr2

  if ( initialized .eq. 0 ) then
    ivec = 0
    ivec(1) = rand_pl_mf
  else
    ivec(1) = rand_pl_mf
    ivec(2) = initialized
    ivec(3) = fib_ctr
    ivec(4) = bd_lastout
    ivec(5) = mersenne_ptr
    ptr1 = 6
    ptr2 = ptr1 + fib_stor_sz - 1
    ivec(ptr1:ptr2) = fib_stor
    ptr1 = ptr2 + 1
    ptr2 = ptr1 + lecuyer_stor_sz - 1
    ivec(ptr1:ptr2) = lecuyer_stor
    ptr1 = ptr2 + 1
    ptr2 = ptr1 + scramble_stor_sz - 1
    ivec(ptr1:ptr2) = scramble_stor
    ptr1 = ptr2 + 1
    ptr2 = ptr1 + mersenne_stor_sz - 1
    ivec(ptr1:ptr2) = mersenne_stor
  end if

end subroutine get_rand_pl_state_main


!=======================================================================
!
!  subroutine set_rand_pl_state_multi(m, ivec)
!  
!  I: ivec   int(rpk) array     contains sufficient information to
!                               restore the state of the generator,
!                               from previous call to get_rand_pl_state
!     m      integer            selects which generator to use
!                               
!  O: only through globals
!  
!  Multiple generator wrapper for set_rand_pl_state_main.
!  
!=======================================================================

subroutine set_rand_pl_state_multi(m, ivec)
  implicit none
  integer(rpk), dimension(rand_pl_state_sz), intent(in) :: ivec
  integer, intent(in) :: m
  if ( curr_ptr .ne. m ) call associate_multi(m)
  call set_rand_pl_state_main(ivec)
end subroutine set_rand_pl_state_multi


!=======================================================================
!
!  subroutine set_rand_pl_state_main(ivec)
!
!  I: ivec   int(rpk) array     contains sufficient information to
!                               restore the state of the generator,
!                               from previous call to get_rand_pl_state
!
!  O: only through globals
!
!=======================================================================

subroutine set_rand_pl_state_main(ivec)

  implicit none

  integer(rpk), dimension(rand_pl_state_sz), intent(in) :: ivec
  integer :: ptr1, ptr2

  if ( ivec(1) .ne. rand_pl_mf ) then
    write(0,*) 'Error (SET_RAND_PL_STATE): rand_pl_mf changed since ' &
               // 'state vector was taken'
    stop
  else
    initialized   = ivec(2)
    fib_ctr       = ivec(3)
    bd_lastout    = ivec(4)
    mersenne_ptr  = ivec(5)
    ptr1 = 6
    ptr2 = ptr1 + fib_stor_sz - 1
    fib_stor      = ivec(ptr1:ptr2)
    ptr1 = ptr2 + 1
    ptr2 = ptr1 + lecuyer_stor_sz - 1
    lecuyer_stor  = ivec(ptr1:ptr2)
    ptr1 = ptr2 + 1
    ptr2 = ptr1 + scramble_stor_sz - 1
    scramble_stor = ivec(ptr1:ptr2)
    ptr1 = ptr2 + 1
    ptr2 = ptr1 + mersenne_stor_sz - 1
    mersenne_stor = ivec(ptr1:ptr2)
  end if

end subroutine set_rand_pl_state_main


end module sderk90