MODULE ran_num
   !!======================================================================
   !!                       *** MODULE ran_num ***
   !! NEMOVAR: Random number routines
   !!======================================================================
   !!----------------------------------------------------------------------
   !! gaustb    : returns a gaussian randum number with mean and std.
   !! gaustb_2d : returns a gaussian randum number with mean and std for a
   !!             given horizontal grid point
   !!----------------------------------------------------------------------
   !! * Modules used
   USE par_kind       ! Kind variables
   USE dom_oce        ! Domain variables
   USE in_out_manager ! I/O stuff

   IMPLICIT NONE

   !! * Routine accessibility

   PRIVATE

   PUBLIC &
      & gaustb, &
      & gaustb_2d

CONTAINS

   FUNCTION gaustb( kseed, pamp, pmean  )
      !!----------------------------------------------------------------------
      !!               ***  ROUTINE gaustb ***
      !!          
      !! ** Purpose : Returns a gaussian randum number with mean and std.
      !!
      !! ** Method  : Generate Gaussian random variables.
      !!              The standard deviation and mean of the variables are 
      !!              specified by the variables pamp and pmean.
      !!
      !! ** Action  : 
      !!
      !! History :
      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
      !!----------------------------------------------------------------------
      !! * Function return
      REAL(wp) :: &
         & gaustb
      !! * Arguments
      INTEGER, INTENT(INOUT) :: &
         & kseed        ! Seed
      REAL(wp), INTENT(IN) :: &
         & pamp, &      ! Amplitude
         & pmean        ! Mean value
      !! * Local declarations
      
      gaustb = pamp * gausva( kseed ) + pmean

   END FUNCTION gaustb


   FUNCTION gausva( kdum )
      !!----------------------------------------------------------------------
      !!               ***  ROUTINE gausva ***
      !!          
      !! ** Purpose : Returns a normally distributed deviate with 0 mean 
      !!              and unit variance using the unifva(kdum) as the 
      !!              source of uniform deviates.
      !!
      !! ** Method  : 
      !!
      !! ** Action  : 
      !!
      !! History :
      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
      !!----------------------------------------------------------------------
      !! * Function return
      REAL(wp) :: &
         & gausva
      !! * Arguments
      INTEGER, INTENT(INOUT) :: &
         & kdum          ! Seed
      !! * Local declarations
      REAL(wp), SAVE :: &
         & gset
      INTEGER, SAVE :: &
         & niset = 0
      REAL(wp) :: &
         & zfac, &
         & zrsq, &
         & zv1,  &
         & zv2

      ! Begin main

      IF ( niset == 0 ) THEN

         zv1   = 2.0_wp * psrandom( kdum ) - 1.0_wp
         zv2   = 2.0_wp * psrandom( kdum ) - 1.0_wp
         zrsq  = zv1**2 + zv2**2
         
         DO WHILE ( ( zrsq >= 1.0_wp ) .OR. ( zrsq == 0.0_wp ) ) 

            zv1   = 2.0_wp * psrandom( kdum ) - 1.0_wp
            zv2   = 2.0_wp * psrandom( kdum ) - 1.0_wp
            zrsq  = zv1**2 + zv2**2

         END DO

         zfac   = SQRT( -2.0_wp * LOG( zrsq ) / zrsq )
         gset   = zv1 * zfac
         gausva = zv2 * zfac
         niset  = 1

      ELSE

         gausva = gset
         niset  = 0

      ENDIF
       
   END FUNCTION gausva

   FUNCTION psrandom( kdum ) 
      !!----------------------------------------------------------------------
      !!               ***  ROUTINE psrandom ***
      !!          
      !! ** Purpose : Pseudo-Random number generator.
      !!
      !! ** Method  : Returns a pseudo-random number from a uniform distribution
      !!              between 0 and 1
      !!              Call with kdum a negative integer to initialize.
      !!              Thereafter, do not alter kdum between successive deviates 
      !!              in sequence.
      !!
      !! ** Action  : 
      !!
      !! History :
      !!        !  10-02  (F. Vigilant)  Original code 
      !!----------------------------------------------------------------------
      !! * Function return
      REAL(wp) ::  &
         & psrandom
      !! * Arguments
      INTEGER, INTENT(INOUT) :: &
         & kdum          ! Seed
      LOGICAL, SAVE :: &
         & llinit = .FALSE.

      ! Initialization
      IF ( .NOT. llinit ) THEN
         
         CALL srand( kdum )
         llinit     = .TRUE.
         
      ENDIF

      CALL rand(psrandom) 

   END FUNCTION psrandom
   

   FUNCTION gaustb_2d( ki, kj, kseed, pamp, pmean  )
      !!----------------------------------------------------------------------
      !!               ***  ROUTINE gaustb_2d ***
      !!          
      !! ** Purpose : Returns a Gaussian randum number with mean and std
      !!              for a given horizontal grid point.
      !!
      !! ** Method  : Generate Gaussian random variables.
      !!              The standard deviation and mean of the variables are 
      !!              specified by the variables pamp and pmean.
      !!
      !! ** Action  : 
      !!
      !! History :
      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
      !!----------------------------------------------------------------------
      !! * Function return
      REAL(wp) :: &
         & gaustb_2d
      !! * Arguments
      INTEGER, INTENT(IN) :: &
         & ki, &        ! Indices in seed array
         & kj
      INTEGER, INTENT(INOUT), DIMENSION(jpi,jpj) :: &
         & kseed        ! Seed
      REAL(wp), INTENT(IN) :: &
         & pamp, &      ! Amplitude
         & pmean        ! Mean value
      !! * Local declarations
      
      gaustb_2d = pamp * gausva_2d( ki, kj, kseed ) + pmean

   END FUNCTION gaustb_2d

   FUNCTION gausva_2d( ki, kj, kdum )
      !!----------------------------------------------------------------------
      !!               ***  ROUTINE gausva_2d ***
      !!          
      !! ** Purpose : Returns a normally distributed deviate with 0 mean 
      !!              and unit variance.
      !!
      !! ** Method  : 
      !!
      !! ** Action  : 
      !!
      !! History :
      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
      !!----------------------------------------------------------------------
      !! * Function return
      REAL(wp) :: &
         & gausva_2d
      !! * Arguments
      INTEGER, INTENT(IN) :: &
         & ki, &         ! Indices in seed array
         & kj
      INTEGER, INTENT(INOUT), DIMENSION(jpi,jpj) :: &
         & kdum          ! Seed
      !! * Local declarations
      REAL(wp), SAVE, DIMENSION(jpi,jpj) :: &
         & gset
      INTEGER, SAVE, DIMENSION(jpi,jpj) :: &
         & niset
      LOGICAL, SAVE :: &
         & llinit = .FALSE.
      REAL(wp) :: &
         & zfac, &
         & zrsq, &
         & zv1,  &
         & zv2

      ! Initialization

      IF ( .NOT. llinit ) THEN

         niset(:,:) = 0
         llinit     = .TRUE.

      ENDIF

      ! Begin main

      IF ( niset(ki,kj) == 0 ) THEN

         zv1   = 2.0_wp * psrandom_2d( ki, kj, kdum ) - 1.0_wp
         zv2   = 2.0_wp * psrandom_2d( ki, kj, kdum ) - 1.0_wp
         zrsq  = zv1**2 + zv2**2

         DO WHILE ( ( zrsq >= 1.0_wp ) .OR. ( zrsq == 0.0_wp ) ) 

            zv1   = 2.0_wp * psrandom_2d( ki, kj, kdum ) - 1.0_wp
            zv2   = 2.0_wp * psrandom_2d( ki, kj, kdum ) - 1.0_wp
            zrsq  = zv1**2 + zv2**2

         END DO

         zfac   = SQRT( -2.0_wp * LOG( zrsq ) / zrsq )
         gset(ki,kj) = zv1 * zfac
         gausva_2d   = zv2 * zfac
         niset(ki,kj) = 1

      ELSE

         gausva_2d = gset(ki,kj)
         niset(ki,kj)  = 0

      ENDIF
       
   END FUNCTION gausva_2d

   FUNCTION psrandom_2d( ki, kj, kdum ) 
      !!----------------------------------------------------------------------
      !!               ***  ROUTINE psrandom_2d ***
      !!          
      !! ** Purpose : Random number generator.
      !!
      !! ** Method  : Returns a pseudo-random number from a uniform distribution
      !!              between 0 and 1
      !!              Call with kdum a negative integer to initialize.
      !!              Thereafter, do not alter kdum between successive deviates 
      !!              in sequence.
      !!
      !! ** Action  : 
      !!
      !! History :
      !!        !  10-02  (F. Vigilant)  Original code
      !!----------------------------------------------------------------------
      !! * Function return
      REAL(wp) ::  &
         & psrandom_2d
      !! * Arguments
      INTEGER, INTENT(IN) :: &
         & ki, &         ! Indices in seed array
         & kj
      INTEGER, INTENT(INOUT), DIMENSION(jpi,jpj) :: &
         & kdum          ! Seed
      LOGICAL, SAVE :: &
         & llinit = .FALSE.

      ! Initialization
      IF ( .NOT. llinit ) THEN
         
         CALL srand( kdum( ki,kj ) )
         llinit     = .TRUE.
         
      ENDIF

      CALL rand(psrandom_2d)  
      
   END FUNCTION psrandom_2d

END MODULE ran_num