!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

module random_mod

  use com_mod, only: maxpart
  ! read: maxpart
  ! write:

  implicit none

  private

  integer, parameter :: ran1_ntab=32
  integer :: ran1_iv(ran1_ntab) = 0, ran1_iy = 0

  integer :: gasdev_iset = 0
  real    :: gasdev_gset

  integer :: ran3_iff = 0
  integer :: ran3_inext,ran3_inextp
  integer :: ran3_ma(55)
  real    :: ran3_initialize(maxpart), ran3_advance(maxpart), ran3_conv(maxpart)

  integer :: idummy_ran3_initialize = -7, idummy_ran3_advance = -7, idummy_ran3_conv = -7
  ! idummy_ran3_initialize, idummy_ran3_advance, idummy_ran3_conv  Initial value for calculation
  !                        of random numbers (see timemanager_computations, convmix_tiedtke,
  !                        and ran3)

  public :: ran1, gasdev, gasdev1, ran3
  public :: ran3_initialize, ran3_advance, ran3_conv
  public :: idummy_ran3_initialize, idummy_ran3_advance, idummy_ran3_conv

contains
  !  Taken from Press et al., Numerical Recipes

  function ran1(idum)
    implicit none
    integer :: idum
    real :: ran1
    integer,parameter :: ia=16807, im=2147483647, iq=127773, ir=2836
    integer,parameter :: ndiv=1+(im-1)/ran1_ntab
    real,parameter :: am=1./real(im), eps=1.2e-7,rnmx=1.-eps
    integer :: j,k
    if (idum.le.0.or.ran1_iy.eq.0) then
      idum=max(-idum,1)
      do j=ran1_ntab+8,1,-1
        k=idum/iq
        idum=ia*(idum-k*iq)-ir*k
        if (idum.lt.0) idum=idum+im
        if (j.le.ran1_ntab) ran1_iv(j)=idum
      enddo
      ran1_iy=ran1_iv(1)
    endif
    k=idum/iq
    idum=ia*(idum-k*iq)-ir*k
    if (idum.lt.0) idum=idum+im
    j=1+ran1_iy/ndiv
    ran1_iy=ran1_iv(j)
    ran1_iv(j)=idum
    ran1=min(am*real(ran1_iy),rnmx)
    return
  end function ran1

  function gasdev(idum)
    implicit none
    integer :: idum
    real :: gasdev, fac, v1, v2, r
    if (gasdev_iset.eq.0) then
1     v1=2.*ran3(idum)-1.
      v2=2.*ran3(idum)-1.
      r=v1**2+v2**2
      if(r.ge.1.0 .or. r.eq.0.0) go to 1
      fac=sqrt(-2.*log(r)/r)
      gasdev_gset=v1*fac
      gasdev=v2*fac
      gasdev_iset=1
    else
      gasdev=gasdev_gset
      gasdev_iset=0
    endif
    return
  end function gasdev

  subroutine gasdev1(idum,random1,random2)
    implicit none
    integer :: idum
    real :: random1, random2, v1, v2, r, fac
1     v1=2.*ran3(idum)-1.
    v2=2.*ran3(idum)-1.
    r=v1**2+v2**2
    if(r.ge.1.0 .or. r.eq.0.0) go to 1
    fac=sqrt(-2.*log(r)/r)
    random1=v1*fac
    random2=v2*fac
    ! Limit the random numbers to lie within the interval -3 and +3
    !**************************************************************
    if (random1.lt.-3.) random1=-3.
    if (random2.lt.-3.) random2=-3.
    if (random1.gt.3.) random1=3.
    if (random2.gt.3.) random2=3.
    return
  end subroutine gasdev1

  function ran3(idum)
    implicit none
    integer :: idum
    real :: ran3
    integer, parameter :: mbig=1000000000, mseed=161803398, mz=0
    real, parameter :: fac=1./real(mbig)
    integer :: i,ii,k
    integer :: mj,mk

    if(idum.lt.0.or.ran3_iff.eq.0)then
      ran3_iff=1
      mj=mseed-iabs(idum)
      mj=mod(mj,mbig)
      ran3_ma(55)=mj
      mk=1
      do i=1,54
        ii=mod(21*i,55)
        ran3_ma(ii)=mk
        mk=mj-mk
        if(mk.lt.mz)mk=mk+mbig
        mj=ran3_ma(ii)
      end do
      do k=1,4
        do i=1,55
          ran3_ma(i)=ran3_ma(i)-ran3_ma(1+mod(i+30,55))
          if(ran3_ma(i).lt.mz)ran3_ma(i)=ran3_ma(i)+mbig
        end do
      end do
      ran3_inext=0
      ran3_inextp=31
      idum=1
    endif
    ran3_inext=ran3_inext+1
    if(ran3_inext.eq.56)ran3_inext=1
    ran3_inextp=ran3_inextp+1
    if(ran3_inextp.eq.56)ran3_inextp=1
    mj=ran3_ma(ran3_inext)-ran3_ma(ran3_inextp)
    if(mj.lt.mz)mj=mj+mbig
    ran3_ma(ran3_inext)=mj
    ran3=real(mj)*fac
    return
  end function ran3
  !  (C) Copr. 1986-92 Numerical Recipes Software us.

end module random_mod