[子程序] fortran程序报错讨论

subroutinernde(ua,n)
 c           random number generator proposed bymarsaglia and zaman
 c           in reportfsu-scri-87-50
 c           modified by f. james, 1988and 1989, to generate a vector
 c           of pseudorandom numbers ua oflength n.
 c           modified by k.k
 c    ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 c     calling sequences for rnde:                                    ++
 c         call rnde(ua, n)         returns a vector ua of n          ++
 c                      32-bit randomfloating point numbers between  ++
 c                      zero and one.                                 ++
 c         call rnd3i(i1)          initializes the generator from one++
 c                      32-bit integeri1
 c   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       implicit none
       integer n
       real*8 ua(n), sina(102)
       logical first
       save  first
       integer ijkl, ijklin
       data first/.true./
 c
       if(first) then
          first=.false.
 c           default initialization. userhas called rnde without rnd3i.
          ijkl = 54217137




c     %%%%%%%%%%%%%%%%%%%%%%%
          call rnd3ix(ijkl)

c    %%%%%%%%%%%%%%%%%%%%%%%%%%%




      endif



c     %%%%%%%%%%%%%%%%%%%%%%%%%

      call rnd3x(ua, n)   

c     %%%%%%%%%%%%%%%%%%%%%%%%%%%




      return
 c         initializing routine for rnde,may be called before
 c         generating pseudorandom numberswith rnde. the input
 c         values should be in theranges:  0<=ijklin<=900 oooooo
 c     **************
       entry rnd3i(ijklin)
 c     *************
       first=.false.
       call rnd3ix(ijklin)
       return
 c     ************
       entry rnd3r(sina)
 c     ************
          first=.false.
          call rnd3rx(sina)
       end
       subroutine  rnd3ix(ijkl)
 c             the standard values inmarsaglia's paper, ijkl=54217137
           implicit none
           integer n
           real*8 ua(n), u(97), uni, s, t,zuni
           real*8 sina(102),sout(102)
           integer  jj, m
           integer ns, ijkl, i97, j97, ij,kl, i, j, k, l, ii, ivec
           real*8 twom24, c, cd, cm
           parameter (ns=24,twom24=2.**(-24))
           save c, cd, cm, i97, j97
 c
           ij = ijkl/30082
           kl = ijkl - 30082*ij
           i = mod(ij/177, 177) + 2
           j = mod(ij, 177) + 2
           k = mod(kl/169, 178) + 1
           l = mod(kl, 169)
            do   ii= 1, 97
               s = 0.
               t = .5
                do   jj= 1, ns
                   m = mod(mod(i*j,179)*k,179)
                   i = j
                   j = k
                   k = m
                   l = mod(53*l+1,169)
                   if(mod(l*m,64) .ge. 32)s = s+t
                   t = 0.5*t
                enddo
               u(ii) = s
            enddo
           c = 362436.*twom24
           cd = 7654321.*twom24
           cm = 16777213.*twom24
           i97 = 97
           j97 = 33
           return
 c     ****************
       entry rnd3x(ua, n)
 c     ****************
        do   ivec= 1, n
           uni = u(i97)-u(j97)
           if(uni .lt. 0.)uni=uni+1.
           u(i97) = uni
           i97 = i97-1
           if(i97 .eq. 0) i97=97
           j97 = j97-1
           if(j97 .eq. 0) j97=97
           c = c - cd
           if(c .lt. 0.) c=c+cm
           uni = uni-c
           if(uni .lt. 0.)uni=uni+1.
           ua(ivec) = uni
 c                 replace exact zeros byuniform distr. *2**-24
           if(uni .eq. 0.) then
               zuni = twom24*u(2)
 c               an exact zero here isvery unlikely, but let's be safe.
               if(zuni .eq. 0.) zuni=twom24*twom24
               ua(ivec) = zuni
           endif
        enddo
       return
 c     ****************** to get currentstatus
       entry rnd3s(sout)
 c     ***********
            do   i=1, 97
               sout(i)=u(i)
            enddo
           sout(98)=c
           sout(99)=cd
           sout(100)=cm
           sout(101)=i97
           sout(102)=j97
           return
 c     ****************  to restore the old status
       entry rnd3rx(sina)
            do   i=1, 97
                u(i)=sina(i)
            enddo
           c=sina(98)
           cd=sina(99)
           cm=sina(100)
           i97=sina(101)
           j97=sina(102)
       end