MODULE GENERATOR
!     =============================================================================
      
      INTEGER, PRIVATE :: DIM
      DOUBLE PRECISION, PRIVATE :: PAR1, PAR2  
      DATA DIM /3/  
!     PAR1 AND PAR2 ARE CONTROL PARAMETERS THAT'S GIVEN BY THE USER.
!     DEFAULT DIMENSION = 3.
      
      CONTAINS
 
      SUBROUTINE ORIG_FIELD (NDIM,XX,X,Y,Z,T,VV,VX,VY,VZ,PARM1,PARM2)
!     ---------------------------------------------------------------------------
!     EVALUATE DESIRED FIELD VARIABLES INSIDE A NDIM-CUBIC BOX OF [0,1].
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER, OPTIONAL, INTENT(IN) :: NDIM
      DOUBLE PRECISION, OPTIONAL, INTENT(IN)  :: XX(:),X,Y,Z,PARM1,PARM2
      DOUBLE PRECISION, OPTIONAL, INTENT(OUT) :: T,VV(:),VX,VY,VZ
      DOUBLE PRECISION PX,PY,PZ,PT,PVX,PVY,PVZ  ! INTERNAL VARIABLES
      
      PX = 0.5D0
      PY = 0.5D0
      PZ = 0.5D0
      IF (PRESENT(XX)) THEN 
      IF (NDIM.GE.1) PX = XX(1)
      IF (NDIM.GE.2) PY = XX(2)
      IF (NDIM.GE.3) PZ = XX(3)
      ELSE
      IF (PRESENT(X)) PX = X
      IF (PRESENT(Y)) PY = Y
      IF (PRESENT(Z)) PZ = Z
      ENDIF      
      
      PVX = 0.D0
      PVY = 0.D0
      PVZ = 0.D0
      IF (PRESENT(NDIM)) DIM=NDIM
      IF (PRESENT(PARM1)) PAR1 = PARM1
      IF (PRESENT(PARM2)) PAR2 = PARM2

!     ------------------------------------------------------------
      PT  = PAR1+PAR2*PX
      PVX = PAR2*PY
!     ------------------------------------------------------------

      IF (PRESENT(T))  T = PT
      IF (PRESENT(VX)) VX = PVX
      IF (PRESENT(VY)) VY = PVY
      IF (PRESENT(VZ)) VZ = PVZ
      IF (PRESENT(VV)) THEN 
      IF (NDIM.GE.1) VV(1) = PVX 
      IF (NDIM.GE.2) VV(2) = PVY
      IF (NDIM.GE.3) VV(3) = PVZ
      ENDIF   
      END SUBROUTINE ORIG_FIELD
!     ---------------------------------------------------------------------------
      
      
      FUNCTION RANDSAFE()
!     ---------------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION RANDSAFE
      SAVE ISEED
      DATA ISEED /1/
 54   RANDSAFE = RAN1(ISEED)
      IF (RANDSAFE.LE.0) GOTO 54
      IF (RANDSAFE.GE.1) GOTO 54
      END FUNCTION RANDSAFE
!     ---------------------------------------------------------------------------

      
      FUNCTION RAN1 (IDUM)
!     ---------------------------------------------------------------------------
      INTEGER IDUM,IA,IM,IQ,IR,NTAB,NDIV
      DOUBLE PRECISION AM,EPS,RNMX,RAN1
      parameter (ia=16807,im=2147483647,am=1.d0/im,iq=127773,ir=2836,ntab=32, & 
      ndiv=1+(im-1)/ntab,eps=1.2d-7,rnmx=1.d0-eps)
      integer j,k,iv(ntab),iy
      save iv,iy
      data iv /ntab*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
        idum=max(-idum,1)
        do 11 j=ntab+8,1,-1
          k=idum/iq
           idum=ia*(idum-k*iq)-ir*k
          if (idum.lt.0) idum=idum+im
          if (j.le.ntab) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/iq
      idum=ia*(idum-k*iq)-ir*k
      if (idum.lt.0) idum=idum+im
      j=1+iy/ndiv
      iy=iv(j)
      iv(j)=idum
      ran1=min(am*iy,rnmx)
      END FUNCTION RAN1
!     ---------------------------------------------------------------------------
      
      END MODULE GENERATOR
!     =============================================================================