c**********************************************************************
c     This block of subroutines contains the additional subroutines
c     and functions needed to run the stochastically forced version
c      of the model.
c********************************************************************
      function urand()
c
c     Return the next pseudo-random deviate from a sequence which is
c     uniformly distributed in the interval [0,1]
c
c     Uses the function ran2, from Press et al, Numerical Recipes,
c     2nd ed., Cambridge Univ. Press, 1992.
c
      implicit none
c
c     Input - none
c
c     Output
      real     urand
c
c     Local
      integer  iseed
      real     ran2
      external ran2
c
c     Common block to make iseed visible to urand_init (and to save
c     it between calls)
      common /urand_seed/ iseed
c
      urand = ran2( iseed )
      return
      end
c*******************************************************************
      subroutine urand_init( urand_seed )
c
c     Initialize random number generator urand with given seed
c
      implicit none
c
c     Input
      integer urand_seed
c
c     Output - none
c
c     Local
      integer iseed
c
c     Common block to communicate with urand
      common /urand_seed/ iseed
c
c     Seet the seed value
      iseed = urand_seed
      return
      end
c******************************************************************
      FUNCTION ran2(idum)
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      REAL ran2,AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 11 j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
      ran2=min(AM*iy,RNMX)
      return
      END