       subroutine RIS(m,nc,indexfd,indexftptm,amu,risamu,iseed)
       implicit real*8 (a-h,o-z)
       real*8 risamu(1:m,1:4,1:9)
       real*8 amu(1:12,1:4,1:9)
       integer nc(4),indexftptm(1:12,1:4,1:9)
       integer indexfd(12)
       integer m,iseed
       real*8 anu
       real*8  ta(4),sig2,tau(1:4,1:9) 
       ta(1)=0.456979d0
       ta(2)=0.428031d0
       ta(3)=0.400165d0
       ta(4)=0.366237d0
       sig2=0.001662d0 
       anu=40.0d0
       do 300 i=1,500
        call Gtau(nc,indexftptm,anu,amu,ta,sig2,tau,iseed)
        if (i .eq. 1) then
        write(*,*) 'i=',i,'tau(1,1)=',tau(1,1)
        endif
        call Gamu1(nc,indexfd,indexftptm,ta,sig2,
     1                  tau,amu,iseed)
 300   continue
c
       icount=0
       do 230 i=1,m
        if ( (i .eq. 1) .or. (i .eq. icount) ) then
         icount=icount+200
        endif
        do 240 j=1,20
         call Gtau(nc,indexftptm,anu,amu,ta,sig2,tau,iseed)
         call Gamu1(nc,indexfd,indexftptm,ta,sig2,
     1                  tau,amu,iseed)
 240    continue
        do 260 i1=1,4
         do 270 j1=1,nc(i1)
          risamu(i,i1,j1)=amu(2,i1,j1)
 270     continue
 260    continue
 230   continue
       return
       end


       subroutine Gtau(nc,indexftptm,anu,amu,ta,sig2,tau,iseed)
       implicit real*8 (a-h,o-z)
       real*8 amu(1:12,1:4,1:9)
       integer nc(4),indexftptm(1:12,1:4,1:9)
       real*8 anu
       real*8  ta(4),sig2,tau(1:4,1:9)
       real*8 r(1),shape,scaler
       integer iseed
       external DRNGAM

       do 200 i=1,4
        if (i .eq. 1) then
         do 210 j=1,nc(i)
          if (indexftptm(2,i,j) .eq. 1) goto 210
          scaler=anu/2.0d0+(amu(2,i,j)-ta(i))**2/(2.0d0*sig2)
          shape=(anu+1.0d0)*0.50d0
          call rnset( iseed )
           call DRNGAM(1,shape,r)
          call rnget( iseed )
          tau(i,j)=r(1)/scaler
 210     continue
        else
         do 220 j=2,nc(i)
          if (indexftptm(2,i,j) .eq. 1) goto 220
          scaler=anu/2.0d0+(amu(2,i,j)-ta(i))**2/(2.0d0*sig2)
          shape=(anu+1.0d0)*0.50d0
          call rnset( iseed )
           call DRNGAM(1,shape,r)
          call rnget( iseed )
          tau(i,j)=r(1)/scaler
 220     continue
        endif
 200   continue
       return
       end


       subroutine Gamu1(nc,indexfd,indexftptm,ta,sig2,
     1                  tau,amu,iseed)
c......generate amu givenz,theta,sqdelta
       implicit real*8 (a-h,o-z)
       real*8 amu(1:12,1:4,1:9),ta(4),sig2,tau(1:4,1:9)
       real*8 aleft,aright,tn
       real*8 amean,sigma
       integer nc(4),indexftptm(1:12,1:4,1:9),indexfd(12)
       integer iseed
       logical la,lb
c......initialize
       la=.false.
       lb=.false.
       do 410 i=1,12
         if (indexfd(i) .eq. 0) goto 410
         lb=.true.
         amean=ta(1)
         sigma=dsqrt(sig2/tau(1,1))
         aleft=amu(i,1,3)
         aleft=max(aleft,amu(i,2,3))
         do 415 j=3,4
           aleft=max(aleft,amu(i,j,2))
 415     continue
         aleft=(aleft-amean)/sigma
         tn=YTUVN(aleft,aright,la,lb,iseed)
         lb=.false.
c........get amu(i,1,1)
         amu(i,1,1)=amean+sigma*tn
         do 420 j=1,4
            do 430 l=2,nc(j)
c.............get two means and sigmas
              if (indexftptm(i,j,l) .eq. 1) goto 430
              amean=ta(j)
              sigma=dsqrt(sig2/tau(j,l))
c.............get the bounds of amu's
c.............begin first row
              if ((j .eq. 1) .and. (l .eq. 3)) then
                aright=(amu(i,1,1)-amean)/sigma
                aleft=(amu(i,2,3)-amean)/sigma
                aleft=max((amu(i,1,6)-amean)/sigma,aleft)
              endif
              if ((j .eq. 1) .and. (l .eq. 6)) then
                aright=(amu(i,1,3)-amean)/sigma
                aleft=(amu(i,2,6)-amean)/sigma
                aleft=max((amu(i,1,7)-amean)/sigma,aleft)
              endif
              if ( (j .eq. 1) .and. (l .lt. nc(1))
     1             .and. (l .gt. 6) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aleft=(amu(i,j,l+1)-amean)/sigma
                aleft=max(aleft,(amu(i,j+1,l)-amean)/sigma)
              endif
              if ( (j .eq. 1) .and. (l .eq. nc(1)) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aleft=(amu(i,j+1,l)-amean)/sigma
              endif
c.............end first row
c.............begin second row
              if ((j .eq. 2) .and. (l .eq. 3)) then
                aright=(amu(i,1,1)-amean)/sigma
                aright=min(aright,(amu(i,1,3)-amean)/sigma)
                aleft=(amu(i,2,4)-amean)/sigma
                aleft=max((amu(i,3,3)-amean)/sigma,aleft)
              endif
              if ((j .eq. 2) .and. (l .eq. 4)) then
                aright=(amu(i,2,3)-amean)/sigma
                aleft=(amu(i,2,5)-amean)/sigma
                aleft=max((amu(i,3,4)-amean)/sigma,aleft)
              endif
              if ((j .eq. 2) .and. (l .eq. 5)) then
                aright=(amu(i,2,4)-amean)/sigma
                aleft=(amu(i,2,6)-amean)/sigma
                aleft=max((amu(i,3,5)-amean)/sigma,aleft)
              endif
              if ( (j .eq. 2) .and. (l .le. nc(3))
     1             .and. (l .gt. 5) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                aleft=(amu(i,j,l+1)-amean)/sigma
                aleft=max(aleft,(amu(i,j+1,l)-amean)/sigma)
              endif
              if ( (j .eq. 2) .and. (l .lt. nc(2))
     1             .and. (l .gt. nc(3)) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                aleft=(amu(i,j,l+1)-amean)/sigma
              endif
              if ( (j .eq. 2) .and. (l .eq. nc(2)) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                la=.true.
              endif
c.............end second row
c.............begin third row
              if ((j .eq. 3) .and. (l .eq. 2)) then
                aright=(amu(i,1,1)-amean)/sigma
                aleft=(amu(i,3,3)-amean)/sigma
                aleft=max((amu(i,4,2)-amean)/sigma,aleft)
              endif
              if ( (j .eq. 3) .and. (l .le. nc(4))
     1             .and. (l .gt. 2) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                aleft=(amu(i,j,l+1)-amean)/sigma
                aleft=max(aleft,(amu(i,j+1,l)-amean)/sigma)
              endif
              if ( (j .eq. 3) .and. (l .lt. nc(3))
     1             .and. (l .gt. nc(4)) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                aleft=(amu(i,j,l+1)-amean)/sigma
              endif
              if ( (j .eq. 3) .and. (l .eq. nc(3)) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                la=.true.
              endif
c.............end third row
c.............begin forth row
              if ((j .eq. 4) .and. (l .eq. 2)) then
                aright=(amu(i,1,1)-amean)/sigma
                aright=min(aright,(amu(i,j-1,2)-amean)/sigma)
                aleft=(amu(i,4,3)-amean)/sigma
              endif
              if ( (j .eq. 4) .and. (l .lt. nc(4))
     1             .and. (l .gt. 2) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                aleft=(amu(i,j,l+1)-amean)/sigma
              endif
              if ( (j .eq. 4) .and. (l .eq. nc(4)) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                la=.true.
              endif
c.............end forth row
c........get amu(i,j,l)
              tn=YTUVN(aleft,aright,la,lb,iseed)
              la=.false.
              lb=.false.
              amu(i,j,l)=amean+sigma*tn
 430        continue
 420     continue
 410   continue
       return
       end


