*  /sd2g/hmaz/doug1/gibbs/exper/expdfap.f
*                copied from exp3.f  -- here I use p(z|kappa)p(data|z,kappa)
*        to draw kappa -- alternative identification scheme.
*  THis is the code that generated the draws for the paper.!!!!!!!!!
*                    12-4-97
*  Running on 3-27-00
*   ***     ******   ***************    *****************    ****   ****   **
      IMPLICIT REAL * 8 (A-H,L,O-Z)
      DOUBLE PRECISION mnlpr
      double precision KAPPA,KVPRI,KSPRI,klg
      INTEGER T,TM,allrdy
      DIMENSION FACVAR(4)
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /USE/ TWO,TOL
      COMMON /USFPAR/ A(500,5,3),B(500,5,3)
      COMMON /PRIORS/ tvpri(3),tspri(3),SIGPRV,
     c   SIGPRI(5,5),mnlpr(3),siglpr(3),vpris(3),spris(3),
     c   KVPRI(3),KSPRI(3)
      COMMON / STORZ/ xmedld(500,3),apintl(500,3)
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /SOME/ SUM
      allrdy = 12500
      NDRAWS = 5000
      NBURIN = 0
*      iseed = 1234579
*      call rnset(iseed)
        T = 351
      TWO = 2.0 d 0
      TOL = 100.0 d 0 * DMACH(4)
      nu1 = dble(t-1)
      tol = 100. d 0 * dmach(4)
      sp1 = nu1 * .5 d 0
      OPEN (UNIT=21,FILE='done2fap.out')
      OPEN (UNIT=12,FILE='gibbsap2fdpb.out')
      OPEN (UNIT=14,FILE='~/doug1/data/kfbonds')
      OPEN (UNIT=18,FILE='BPARSap2fdpb.out')
      OPEN (UNIT=20,FILE='ALLZSap2fdpb.out')
      OPEN (UNIT=29,FILE='BADZSap2fdpb.out')
      M = 5
      k = 2
* M = number of bonds 
* K = number of factors
* T = length of time series
      nu2 = dble(t-m)
      sp2 = .5 d 0 * nu2
      NITERS = 0
       ktvgt = 0
       nsk = 0
      mm = m - 1
      mp = m + 1
      TM = T + M
      s = 0.0 d 0
*  READ in the data.
      do 2 i = 1,T
 2    READ (14,1466) idate,r(i,1),smt(i,1),r(i,2),smt(i,2),r(i,3),
     . smt(i,3),r(i,4),smt(i,4),r(i,5),smt(i,5)
 1466 format (2x,i6,2x,f6.4,1x,f7.4,2x,f6.4,1x,f7.4,2x,f6.4,1x,
     .f8.4,2x,f6.4,1x,f8.4,2x,f6.4,1x,f9.4)
*     read (14,1466) idate,r(i,1),r(i,2),r(i,3),r(i,4),r(i,5)
*1466 format (i4,5(1x,f9.7) )
*     smt(i,1) = 13.0 / 52.14
*     smt(i,2) = 26.0 / 52.14
*     smt(i,3) = 5.0 d 0
*     smt(i,4) = 15.0 d 0
* 2   smt(i,5) = 25.0 d 0
* * * * * ** * * **** * **** ***** *********    *   *    *   * * ************ 
       do 2201 ju = 1,t
       do 2201 jb = 1,M
 2201  smt(ju,jb) = smt(ju,jb) / 52.14 d 0
      dt = dble(t)
*      if (allrdy.ne.0) then
           read (21,1866) ngdr,(theta(j),j=1,2),(sigf(jg), jg= 1,2),
     r  (kappa(jk), jk=1,2), (LAMBDA(JL),jl=1,2)
*       theta(1) = theta(1) + theta(2)
        DO 2098 im = 1,M
 2098   READ (21,2066) (SIGMA(im,mj), mj = 1,M)
       CALL DLINRG(M,SIGMA,5,SIGINV,5)
        do 2099 iu = 0,T
        READ (21,2267) itt,(z(iu,nfa), nfa=1,2)
*         z(iu,1) = z(iu,1) + z(iu,2)
*       print 678, iu,(z(iu,nfb), nfb = 1,k)
 678  format (i5,3(f12.9) )
*           z(iu,1) = z(iu,1) + .6 * z(iu,3)
*           z(iu,2) = z(iu,2) + .4 * z(iu,3)
           do 2096 ifa = 1,K
 2096       xmedld(iu,ifa) = z(iu,ifa) * .9998 d 0
 2099   CONTINUE
*         end if
* * * * * ** * * **** * **** ***** *********    *   *    *   * * ************
*   ***     ******   ***************    *****************    ****   ****   **
*  Now we start the iterations
*   ***     ******   ***************    *****************    ****   ****   **
        do 2000 ngdr = 1,NDRAWS+nburin
*   ***     ******   ***************    *****************    ****   ****   **
*  draw Z(t) conditional on Everything else:
      DO 101, kf = 1,K
        FACVAR(kf) = sigf(kf) * sigf(kf)
        gft = (kappa(kf) + LAMBDA(KF) ) * (kappa(kf) + LAMBDA(KF) )
        GAMMA(kf) = DSQRT(gft + TWO * FACVAR(KF))
        klg = kappa(kf) + LAMBDA(KF) + GAMMA(KF)
        aexp = (two * kappa(kf) * theta(kf) ) / facvar(kf)
        tg = two * gamma(kf)
      DO 101 it = 1,T
      DO 101 ib = 1,M
       expo1 = gamma(kf)*smt(it,ib)
       expo2 = 0.5d 0 * klg * smt(it,ib)
        if (expo1.lt.600.0 d0.and.expo2.lt.600.0 d0) then
         ant = tg * dexp(expo2)
         adt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two * ( dexp(expo1) - 1.0 d0 )
         bdt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        else if (expo1.ge.600.0 d0.and.expo2.ge.600.0 d0) then
         ant = tg
         adt = klg * dexp(expo1 - expo2)
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two
         bdt = klg
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        else if (expo1.lt.600.0 d0.and.expo2.ge.600.0 d0) then
         ant = tg
         adt = klg *  dexp(expo1 - expo2)
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two * ( dexp(expo1) - 1.0 d0 )
         bdt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        else if (expo1.ge.600.0 d0.and.expo2.lt.600.0 d0) then
         ant = tg * dexp(expo2 - expo1)
         adt = klg
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two
         bdt = klg
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        end if
  101 CONTINUE
        print 8421
 8421  format (' calling augz. ')
        CALL AUGZ
*         print 2658, Z(0,1),Z(0,2),Z(T+1,1),Z(T+1,2)
 2658 format (' z(0): ',2f12.6,' z(T+1): ',2F12.6)
        CALL GETZ
        print 6001
        CALL DRWSIG
        DO 6549 im = 1,M
 6549   PRINT 1655, (SIGMA(im,jm), jm = 1,M)
        print 6003
        CALL GETTH
        print 1656, (theta(jl),jl=1,K)
        print 6002
        CALL GETK
        PRINT 1657, (kappa(jl),jl=1,K)
        print 6004
        CALL GETL
        print 1658, (lambda(jl),jl= 1,K)
        print 6005
        CALL GETSF
         print 1659, (sigf(jl),jl=1,K)
 6001 format (' getz fini ')
 6002 format (' getth fini ')
 6003 format (' drwsig fini ')
 6004 format (' getk fini ')
 6005 format (' getl fini ')
 1654 format (i3,1x,f17.9)
 1655 format (5(1x,f19.11))
 1656 format (' THETA: ',3f19.9)
 1657 format (' KAPPA: ',3f19.9)
 1659 format (' sig fac: ',3f19.9)
 1658 format (' LAMBDA: ',3f19.9)
 2651 format (' KAPPA: ',3f15.9)
*       if (ngdr.eq.nburin) then
*           isa = 1
*           go to 2000
*        end if
*        if (isa.eq.1) then
          ngkdr = ngdr + allrdy
          WRITE (18,1866) ngkdr,(theta(j),j=1,k),(sigf(jg), jg= 1,k),
     w  (kappa(jk), jk=1,K), (LAMBDA(JL),jl=1,K)
          do 3999 im = 1,M
 3999 WRITE (12,2066) (SIGMA(im,mj), mj = 1,M)
          write (12,2067) 
 2067 format (/)
*          END IF
       print 2656, ngdr
*       write (26,2666) iseed
 2666 format (i44)
 2656 format (' completed ',i5,' draws.')
*     IF (ngdr.eq.ndraws) then
         do 2001 it = 0,T+1
 2001 write (20,2267) it,(z(it,nfa), nfa=1,k)
*     write (20,2056) 
 2056 format (/)
*      END IF
 2066 format (6(1x,f13.10) )
 2000 CONTINUE
 1866 format (i5,12(1x,f11.8) )
 1999 format (12(1x,f11.8) )
 2267 format (i3,3(1x,f12.9) )
        STOP
        END
        SUBROUTINE GETZ
        IMPLICIT REAL * 8 (A-H,L,O-Z)
        DOUBLE PRECISION mnlpr,kappa,kvpri,kspri
        DIMENSION TABLE1(500),table2(500),udifg(100),XV(0:100000),
     d pval(0:100000)
        INTEGER T
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /PRIORS/ tvpri(3),tspri(3),SIGPRV,
     c   SIGPRI(5,5),mnlpr(3),siglpr(3),vpris(3),spris(3),
     c   KVPRI(3),KSPRI(3)
      COMMON /USE/ TWO,TOL
      COMMON / STORZ/ xmedld(500,3),apderl(500,3)
      COMMON /USFPAR/ A(500,5,3),B(500,5,3)
      COMMON /WHERE/ kneed,ineed,ndraw
      COMMON /SOME/ SUM
      COMMON /LOWBOU/ zlow
      COMMON /ZLN/ ZLRGNO,zln2,zln3,nc
      EXTERNAL DRAWZ,pdfz
        errabs = 0.0 d 0
        zero = 0.0 d 0
        errel = .1 d -1
       do 1 ineed = 1,T
*      print 2565, ineed
       do 1 kneed = 1,K
       nc = 0
       iw = 0
       gnn = pdfz(z(ineed,kneed) )
        difm = 1.0 d 0
        iaborh = 0
        iabort = 0
 5000   ibk = 0
****
        if (ndraw.le.1) then
             pvmax = -.1 d -15
             tpv = -.1 d -15
            xg = z(ineed,kneed) * 0.35 d 0
             diffo = xg / 333.333 d 0
            do 123 ig = 1,10000
            xg = xg + diffo
            xv(ig)  = pdfz(xg)
              if (ig.eq.1) then
                 if (xv(ig).gt.9 d 25) then
                    nc = 0
                    xg = xg * .5 d 0
                    go to 123
                 end if
              end if
*           print 7685, xg,xv(ig)
 7685  format (1x,f7.6,2x,f33.6)
            if (xv(ig).gt.tpv) then
               oldmax = pvmax
               pvmax = xv(ig)
               jmax = ig
               xmxa = pvmax * 0.0005 d 0
               tpv = pvmax - .1 d -6
             else
               if (xv(ig).lt.xmxa) go to 1011
             END IF
 123        CONTINUE
 1011  do 1012 iu = 1,10000
       if (xv(iu).ge.xmxa) then
             trystr = z(ineed,kneed) * 0.35 d 0 + dfloat(iu) * diffo
             t21 = xv(iu)
             trylas = z(ineed,kneed) * 0.35 d 0 + dfloat(ig) * diffo
             go to 33
        end if
 1012 CONTINUE
        END IF
****
        eval1 = pdfz(z(ineed,kneed) )
        eval2 = pdfz(xmedld(ineed,kneed) )
*       print 5643, z(ineed,kneed),gnn,xmedld(ineed,kneed),eval2
 5643 format (' start: ',f6.5,1x,f6.3,2xf6.5,1x,f44.3)
        DIF = DABS(z(ineed,kneed) - xmedld(ineed,kneed) ) * 0.25 d 0
*         if (iabort.ge.1.or.iaborh.ge.1) then
*       print 6890, z(ineed,kneed),xmedld(ineed,kneed),eval1,eval2
*           end if
 6890 format (' 6890: ',4f16.9) 
            if (z(ineed,kneed).gt.xmedld(ineed,kneed) ) then
                 ibiger = 1
            else 
                 ibiger = 0
            end if
        if (eval1.ge.eval2) then
               xkmax = eval1
                tryout = xmedld(ineed,kneed)
             if (ibiger.eq.1) then
 5467           tryout = tryout + dif
                 eval3 = pdfz(tryout)
*            print 5675, tryout,eval3,xkmax
 5675 format (' tryout: ',f9.7,' pdf: ',f37.5,' xkmax: ',f9.7)
                 if (eval3.gt.xkmax) then
                    xkmax = eval3
                    go to 5467
                 end if
                 zstart = tryout - dif
             else 
 5468            tryout = tryout - dif
                 if (tryout.le.zero) then 
                       tryout = tryout + dif
                       dif = dif * .4 d 0
                       go to 5468
                  end if
                  eval3 = pdfz(tryout)
*            print 5675, tryout,eval3,xkmax
                  if (eval3.gt.xkmax) then
                        xkmax = eval3
                        go to 5468
                 end if
                  zstart = tryout + dif
             end if
         else
              xkmax = eval2
            tryout = z(ineed,kneed)
             if (ibiger.eq.1) then 
 5469           tryout = tryout - dif 
                 if (tryout.le.zero) then 
                       tryout = tryout + dif
                       dif = dif * .4 d 0
                       go to 5469
                 end if
                 eval3 = pdfz(tryout) 
*            print 5675, tryout,eval3,xkmax
                 if (eval3.gt.xkmax) then
                     xkmax = eval3
                     go to 5469 
                 end if
                 zstart = tryout + dif
             else 
 5470            tryout = tryout + dif
                  eval3 = pdfz(tryout)
*            print 5675, tryout,eval3,xkmax
                  if (eval3.gt.xkmax) go to 5470 
                  zstart = tryout - dif
             end if
         end if
 7000         dif = .0008 d 0 * difm
             if (z(ineed,kneed).lt.0.008 d 0) dif = dif * .6 d 0
              if (iabort.ge.1.or.iaborh.ge.1) then
                  denom = max0(iabort,iaborh)
                 dif = dif / (denom * 10.0 d 0)
               end if
              trystr = zstart - .00002 d 0
             if (trystr.le.zero) then
                     trystr = .1 d -10 + dif
                end if
              trylas = zstart + .00002 d 0
 4555         trystr = trystr - dif
                 if (trystr.le.zero) then
                     trystr = trystr + dif
                      dif = dif * 0.4 d 0
                      go to 4555
                 end if
              ev = pdfz(trystr)
*            print 5675, trystr,ev,xkmax
              t21 = ev
                if (iabort.ge.1.or.iaborh.ge.1) then
                   print 4546, trystr,ev
 4546 format (' at ',f13.9,' kernel = ',f22.11)
                 end if
              relrat = ev / xkmax
              if (trystr.le..1 d -5) go to 4556
              if (relrat.gt..0005 d 0) go to 4555
 4556         trylas = trylas + dif
              ev = pdfz(trylas)  
*            print 5675, trystr,ev,xkmax
              relrat = ev / xkmax
                if (iabort.ge.1.or.iaborh.ge.1) then 
                   print 4546, trylas,ev
                 end if
              if (relrat.gt..0005 d 0) go to 4556
              trystr = dmax1(trystr,0.0 d 0)
 2566 format (' kneed = ',i1)
 2565 format (' ineed = ',i4)
*       print 9787, ineed,kneed,trystr,trylas
 9787 format (2(1x,i4),' trystr: ',f16.9,' trylas: ',f16.9)
  33   iptsig = 50
       iepts = 78
       difg = (trylas - trystr) / dfloat(iptsig + iepts)
 8000  do 8907 itp = 1,iptsig
 8907 udifg(itp) = difg
       table1(1) = trystr
       table2(1) = t21
*       print 6766, ineed,kneed,table1(1),table2(1)
          udifg(2) = 20.0 d 0 * difg
          udifg(3) = 10.0 d 0 * difg
          udifg(4) = 6.0 d 0 * difg
          udifg(5) = 3.0 d 0 * difg
          udifg(49) = udifg(2)
          udifg(48) = udifg(3)
          udifg(47) = udifg(4)
          udifg(46) = udifg(5)
       do 76 ig = 2,iptsig
       table1(ig) = table1(ig-1) + udifg(ig)
       table2(ig) = pdfz(table1(ig))
*      print 6765, ig,table1(ig),table2(ig)
 6765  format (1x,i4,2x,f10.8,2x,f33.4)
 6766  format (1x,i4,2x,i2,1x,f10.8,2x,f33.4)
  76   CONTINUE
       prob = drnunf()
        oldz = z(ineed,kneed)

***** 
        do 7772, iv = 1, iptsig -1
         if(table1(iv+1).le.table1(iv)) then
           iw = 1
           go to 8888
         end if
 7772  CONTINUE
        z(ineed,kneed) = DGCIN(prob,2,iptsig,table1,table2)
        xmedld(ineed,kneed) = DGCIN(0.5 d 0,2,iptsig,table1,table2)
         IF (Z(INEED,KNEED).GT.0.0 d 0.and.Z(INEED,KNEED).LT.1.0 d 0) 
     i              GO TO 1
*            print 7659, ineed,kneed,z(ineed,kneed)
        zcomp = 0.75 d 0 * z(ineed-1,kneed) + 0.75 d 0 * 
     e          z(ineed+1,kneed)
        ozc = oldz * 1.4 d 0
         if (z(ineed,kneed).gt.zcomp.or.z(ineed,kneed).gt.ozc)
     i      then
             dii = 0.75 d -3
             go to 8002
            end if
        zcomp = 0.2 d 0 * z(ineed-1,kneed) + 0.2 d 0 *
     e          z(ineed+1,kneed)
         zacom = .5 d 0 * z(ineed-1,kneed)
        ozc = oldz * 0.5 d 0   
         if (z(ineed,kneed).lt.zcomp.or.z(ineed,kneed).lt.ozc.or.
     i       z(ineed,kneed).lt.zacom) then
             dii = 0.5 d -4 
             go to 8002 
            end if
         GO TO 1
 8002        print 7659, ineed,kneed,z(ineed,kneed)
 7659 format (' **** z(',i3,',',i1,') ' ,f12.8)
             beyrec = 0.09 d 0
             nts = 0
 8888        mu0 = 0.80 d 0
 8889         umin = .0005 d 0
             xv(0) = mu0 * z(ineed,kneed)
             pes = pdfz(xv(0))
              if (pes.gt.umin) then
                      mu0 = mu0 * .8 d 0
                      go to 8889
              end if
*      print 1088, xv(0),pval(0)
             pvmax = -.1 d -8
            DO 1210 jt = 1,100000
          xv(jt) = xv(jt-1) + dii
            if (xv(jt).gt.beyrec) then
            nts = nts + 1
             dii = .1 d -5 / dfloat(nts)
         ZLRGNO = zlrgno * 2.0 d 0
         zln2 = zln2 * 2.0 d 0 
         zln3 = zln3 * 2.0 d 0
            if (nts.gt.9) then
             xv(0) = 0.0 d 0
             do 1222 jp = 1,10000
          xv(jp) = xv(jp-1) + dii
         pval(jp) = pdfz(xv(jp))
           if (iw.eq.1) then
          write (29,1088) xv(jp),pval(jp)
           end if
 1222    CONTINUE
                  STOP 9823
            end if
         go to 8888
             end if
         pval(jt) = pdfz(xv(jt))
*      print 1088, xv(jt),pval(jt)
 1088 format (f18.16,1x,f18.7)
            tpv = pvmax - .1 d -7
            if (pval(jt).gt.tpv) then
               oldmax = pvmax
               pvmax = pval(jt)
               jmax = jt
               xmxa = pvmax * 0.0005 d 0
 765  format (2(1x,f52.43) )
             else
*            print 6758, xmxa
 6758 format (' xmxa: ',f17.9)
                if (pval(jt).gt.0.0 d 0.and.pval(jt).lt.9 d 9)
     i                       dii = 1.1 d 0 * dii
               if (pval(jt).lt.xmxa) go to 1211
             end if
 1210  CONTINUE
 1211  trylas = xv(jt)
       table1(iptsig) = trylas
       table2(iptsig) = pval(jt)
         DO 102 iu = 1,jt
            if (pval(iu).gt.xmxa) go to 1021
  102    CONTINUE
 1021  table1(1) = xv(iu-1)
       table2(1) = pval(iu-1)
 4563 format (1x,f9.5,2x,f39.37)
       difg = (table1(iptsig) - table1(1)) / dfloat(iptsig + iepts)
          udifg(2) = 20.0 d 0 * difg
          udifg(3) = 10.0 d 0 * difg
          udifg(4) = 6.0 d 0 * difg
          udifg(5) = 3.0 d 0 * difg
          udifg(49) = udifg(2)
          udifg(48) = udifg(3)
          udifg(47) = udifg(4)
          udifg(46) = udifg(5)
       do 276 ig = 2,iptsig
       table1(ig) = table1(ig-1) + udifg(ig)
       table2(ig) = pdfz(table1(ig))
 276   CONTINUE

***** 
        do 7773, iv = 1, iptsig - 1
         if(table1(iv+1).le.table1(iv)) then
           write(50, 7990)
 7990 format (' after manual adj. still bad table')
           write(50, 7991) ineed,kneed
 7991 format (' at pt: ',2i4)
           WRITE (50,1866) (theta(j),j=1,k),(sigf(jg), jg= 1,k),
     .  (kappa(jk), jk=1,K), (LAMBDA(JL),jl=1,K)
 1866 format (12(f12.8) )
      WRITE (50,2066) (SIGMA(im,mj), mj = 1,M)
 2066 format (5(f12.9) )
         do 2002 it = 0,T+1
 2002 write (50,2267) it,(z(it,nfa), nfa=1,k)
 2267 format (i3,3f12.8)
        do 7658 ipt = 1,iptsig
 7658 write (50,5656) ipt,table1(ipt),table2(ipt)
 5656 format(i3,f12.8,2x,f19.9)
          xv(0) = 0.0 d 0
          dii = 0.1 d -5
         do 2176 ip = 1,100000
         xv(ip) = xv(ip-1) + dii
         ph = pdfz(xv(ip))
 2176     write (50,766) xv(ip),ph       
  766 format (f10.8,f26.12)
         STOP 6646
           end if
 7773   continue
        z(ineed,kneed) = DGCIN(prob,2,iptsig,table1,table2)
 1      CONTINUE
        RETURN
        END
      DOUBLE PRECISION FUNCTION PDFZ(x)
*  subroutine to evaluate the kernel of z(ineed,kneed)
*   ineed is the point in time, kneed is the number of factor.
*    (For the 1 factor model, kneed is always 1.)
        IMPLICIT REAL * 8 (A-H,L,O-Z)
        DOUBLE PRECISION mnlpr,kappa,BSI(1000)
        DOUBLE PRECISION kvpri,kspri
        INTEGER T
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /PRIORS/ tvpri(3),tspri(3),SIGPRV,
     c   SIGPRI(5,5),mnlpr(3),siglpr(3),vpris(3),spris(3),
     c   KVPRI(3),KSPRI(3)
      COMMON /USE/ TWO,TOL
      COMMON /USFPAR/ A(500,5,3),B(500,5,3)
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /SOME/ SUM
      COMMON /ZLN/ ZLRGNO,zln2,zln3,nc
      DIMENSION EP(5),t1(5)
         nc = nc + 1
         pi = 3.141592653589793238462643 d 0
        IP = 0
      DO 2, nb = 1,M
       pv = A(ineed,nb,kneed)-B(ineed,nb,kneed)*x
      DO 1, nf = 1,K
          if (nf.eq.kneed) go to 1
         pv = pv + A(ineed,nb,nf) - B(ineed,nb,nf) * Z(ineed,nf)
  1   CONTINUE
      ep(nb) = r(ineed,nb) + pv
  2   CONTINUE
*      do 2626 ib = 1,M
*2626  print 6262, x,r(ineed,ib),ep(ib)
 6262 format (3(1x,f12.4))
         s = 0.0 d 0
      do 22, j = 1,m 
  22     t1(j) = 0.0 d 0
          do 31, ij = 1,M
          do 31, j = 1,M
 31       t1(ij) = t1(ij) + ep(j) *SIGINV(j,ij) 
        do 32, j = 1,m
 32       s = s + t1(j) * ep(j)
       if (nc.eq.1) zlrgno = .5 d 0 * s
      ft = dexp(-.5 d 0 * s + zlrgno)
*  ft is the first term in the conditional density for z, i.e., f(y(t)|z(t)
*  as per the likelihood.
* now we get the second term: f(z(t) | z(t-1)
* Note that here we are working only with 1 of the {z} at a time, whereas in 
*   the preceding portion, we conditioned on the other Z(t) (K > 1)
*  Here I'm using the CIR variables on p. 392.
*        print 7867, ineed,kneed,z(ineed-1,kneed),kappa(kneed),
*    p theta(kneed),sigf(kneed)
 7867 format (1x,i3,1x,i3,4(1x,f12.6) )
        vf = sigf(kneed) * sigf(kneed)
        dt = 1.0 d 0 - DEXP(-KAPPA(KNEED)*(1.0 d 0/52.14 d 0) )
        CIRC = (two * kappa(kneed)) / (vf * dt)
        CIRU = CIRC * z(ineed-1,kneed) * DEXP(-KAPPA(kneed) *
     e         (1.0 d 0/52.14 d 0) )
        CIRV = CIRC * x
        CIRQ = ( (two*kappa(kneed)*theta(kneed))/vf ) - 1.0 d 0
       ACIRQ = dabs(CIRQ)
      bbx = (CIRV / CIRU) ** (.5 d 0 * CIRQ)
      bx = two * ( (CIRU * CIRV)**.5 d0)
       if (ineed.ge.136) then
*         print 4563, bx
           end if
 4563 format (' our bx = ',f33.20)
       if (bx.gt.1450.0 d 0) then
           zec = (8.0 d 0 * bx) * (8.0 d 0 * bx) * (8.0 d 0 * bx)
           BEST = 1.0 d 0 / dsqrt(two*pi*bx)
           bxmu = 4.0 d 0 * acirq*acirq
           best = best * (1.0 d 0 - (bxmu-1.0 d 0)/(8.0 d 0 * bx) +
     e ((bxmu-1.0 d 0)*(bxmu-9.0 d 0))/(two*(8.0 d 0*bx)*(8.0 d 0*bx))
     e - ((bxmu-1.0 d 0)*(bxmu-9.0 d 0)*(bxmu-25.0 d 0))/
     e  (6.0 d 0 *zec))
         if (best.lt.0.0 d 0) then
             best = 1.0 d -6
*            print 2398, best
          end if
 2398 format (' neg best approx used: ',f22.12)
      else
      if (ACIRQ.lt.1.0 d 0) then
         iu = 1
         go to 313
       end if
*     do 310 i = 2,200
*        if (dfloat(i).gt.ACIRQ.and.dfloat(i-1).lt.ACIRQ ) go to 312
*310   CONTINUE
*     Print 3166, ACIRQ
 3166 format (' problem - ACIRQ = ',f12.5)
*  Note x must be positive, and xnu must fall between 0 and 1.  To allow
*  for nu values greater than 1, we have to take the integer amount of 
*  draws.  The last one will be our boy.
* 312  iu = i-1
*      ACIRQ = ACIRQ - dfloat(iu)
       diu = dint(ACIRQ)
       iu = int(diu) + 1
       ACIRQ = DMOD(ACIRQ,diu)
 3136 format (' bx = ',f15.4)
 313   CALL DBSIES(ACIRQ,bx,iu,bsi)
* p. 10 in the special functions manual.
      BEST = BSI(iu)
        END IF
 4565 bby = CIRC * DEXP(-CIRU-CIRV+bx)
      if (nc.eq.1) zln2 = 0.2 d 0 / (bby * bbx * BEST)
      st = (bby * bbx * BEST) * zln2

* NOW get the 3rd term: p(z(t+1) | z(t) )
        CIRU = CIRC * x * DEXP(-KAPPA(kneed) *(1.0 d 0/52.14 d 0) )
        CIRV = CIRC * z(ineed+1,kneed)
      bbx = (CIRV / CIRU) ** (.5 d 0 * CIRQ)
      bx = two * ( (CIRU * CIRV)**.5 d0)
           
       if (bx.gt.1450.0 d 0) then
           zec = (8.0 d 0 * bx) * (8.0 d 0 * bx) * (8.0 d 0 * bx)
           BEST = 1.0 d 0 / dsqrt(two*pi*bx)
           bxmu = 4.0 d 0 * acirq*acirq
           lasn = (bxmu - 1.0 d0) * (bxmu - 9.0 d0) * (bxmu - 25.0 d0)
     * * (bxmu - 49.0 d0)
           lasd = 24.0 d0 * zec * 8.0 d0 * bx
           lasnd = lasn / lasd
*** problemo
           best = best * (1.0 d 0 - (bxmu-1.0 d 0)/(8.0 d 0 * bx) +
     e ((bxmu-1.0 d 0)*(bxmu-9.0 d 0))/(two*(8.0 d 0*bx)*(8.0 d 0*bx))
     e - ((bxmu-1.0 d 0)*(bxmu-9.0 d 0)*(bxmu-25.0 d 0))/
     e  (6.0 d 0 *zec) + lasnd)
         if (best.lt.0.0 d 0) then
              best = 1.0 d -6
*            print 2398, best
          end if
*** problemo
      else
          if (ACIRQ.lt.1.0 d 0) then
             iu = 1
             go to 2313
           end if
*         do 2310 i = 2,200
*         if (dfloat(i).gt.ACIRQ.and.dfloat(i-1).lt.ACIRQ ) go to 2312
*2310     CONTINUE
*         Print 3166, ACIRQ
*  Note x must be positive, and xnu must fall between 0 and 1.  To allow
*  for nu values greater than 1, we have to take the integer amount of
*  draws.  The last one will be our boy.
*2312       iu = i-1
*           ACIRQ = ACIRQ - dfloat(iu)
       diu = dint(ACIRQ)
       iu = int(diu) + 1
       ACIRQ = DMOD(ACIRQ,diu)
 2313       CALL DBSIES(ACIRQ,bx,iu,bsi)
* p. 10 in the special functions manual.
            BEST = BSI(iu)

            if (ineed.ge.136) then
             do jju = 1, iu
*            print 4183, jju, bsi(jju)
*            print 4184, best
             end do   
 4183       format('bsi ',i2,' = ',f9.7)
 4184       format('best ',' = ',f9.7)
            end if
      end if

      bby = CIRC * DEXP(-CIRU-CIRV+bx)
       if (nc.eq.1) zln3 = 0.9 d 0 / (bby * bbx * BEST)
      tt = (bby * bbx * BEST) * zln3

       if (ineed.ge.136) then
*          print 3265, bby,bbx,BEST
       end if
 3265 format (' ttt: ',3f19.12)


      pdfz = (ft * st * tt)
*     pdfz = (ft * st * tt) / sum
        IF (pdfz.lt.0.0 d 0) then 
      PRINT 2366, x,ft,st,tt,sum
        STOP 5675
        END IF
 2366 format (' negat pdf: ',f8.6,3(1x,f17.7),1x,f12.3)
      RETURN
      END

      SUBROUTINE GETTH
      IMPLICIT REAL * 8 (A-H,L,O-Z)
      DOUBLE PRECISION mnlpr,kappa,kvpri,kspri
        DIMENSION TABLE1(500),TABLE2(500),xv(0:10000),pval(10000),
     d   udifg(900)
      INTEGER T
      COMMON /DATA/ R(500,5),Z(0:501,3),smt(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     .   sigf(3),SIGINV(5,5)
      COMMON /PRIORS/ tvpri(3),tspri(3),sigprv,
     .   WW(5,5),mnlpr(3),siglpr(3),vpris(3),spris(3),
     .   KVPRI(3),KSPRI(3)
      COMMON /use/ two,tol
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /TCONS/ xtkons,zlnt,ntimth
      EXTERNAL pdfth
       diffi = 0.07 d -2
       iptsig = 90
        errabs = 0.0 d 0
        zero = 0.0 d 0
        errel = .1 d -1
       do 1 kneed = 1,K
          diffo = theta(kneed) * .03 d -1
*         if (kneed.eq.2) diffo = diffo * 0.33 d 0
          pvmax = -1.0 d -6
 2347      ntimth = 0

*        this is to find th constant
          garb = pdfth(theta(kneed))
          xv(0) = 0.92 d 0 * theta(kneed)
*         if (kneed.eq.1) xv(0) = .9 * theta(kneed)
          JNEED = 10000
 4565      if (xv(0).lt.zero) xv(0) = zero
          pvmax = -1.0 d -6
         DO 1010 jt = 1,JNEED
          xv(jt) = xv(jt-1) + diffo
          pval(jt) = pdfth(xv(jt))
*         print 6758, xv(jt),pval(jt)
 6758 format (1x,f9.6,1x,f20.10)
            tpv = pvmax - .1 d -4
            if (pval(jt).gt.tpv) then
               oldmax = pvmax
               pvmax = pval(jt)
               jmax = jt
               xmxa = pvmax * 0.0005 d 0
 765  format (2(1x,f52.43) )
             else
               if (pval(jt).lt.xmxa) go to 1011
             END IF
 1010  CONTINUE 
 1011  trylas = xv(jt)
*      print 6378
 6378  format(' I got to this point ')
       table1(iptsig) = trylas
       table2(iptsig) = pval(jt)
         DO 102 iu = 1,jt
            if (pval(iu).gt.xmxa) then
               if (iu.eq.1) then
               print 7865
 7865 format (' starting point was too large.')
                  xv(0) = .8 d 0 * xv(0)
                  GO TO 4565
               end if
              go to 1021
             end if
  102    CONTINUE
*      print 6378
 1021  trystr = xv(iu-1)
       table2(1) = pval(iu-1)
       iptsig = 90
       iepts = 78
       difg = (trylas - trystr) / dfloat(iptsig + iepts)
        if (difg.le.zero) then
           difg = dabs(difg)
*           print 1026, trylas,trystr,difg
        end if
 1026 format (' negative difg. trylas: ',f12.8,' trystr: ',f12.7,' difg:
     f   ',f9.6)    
 8000  do 8907 itp = 1,iptsig
 8907 udifg(itp) = difg
       table1(1) = trystr
          udifg(2) = 20.0 d 0 * difg
          udifg(3) = 10.0 d 0 * difg
          udifg(4) = 6.0 d 0 * difg
          udifg(5) = 3.0 d 0 * difg
          udifg(iptsig-1) = udifg(2)
          udifg(iptsig-2) = udifg(3)
          udifg(iptsig-3) = udifg(4)
          udifg(iptsig-4) = udifg(5)
       do 76 ig = 2,iptsig
       table1(ig) = table1(ig-1) + udifg(ig)
       table2(ig) = pdfth(table1(ig))
*       PRINT 7676, TABLE1(IG),TABLE2(IG)
 7676 format (2(1x,f22.6) )
  76   CONTINUE
       prob = drnunf()
 
*****
        do 7772, iv = 1, iptsig - 1
         if(table1(iv+1).le.table1(iv)) then
           write(50, 7990)
 7990   format(' problem in theta ')
           write(50, 7991) kneed
 7991   format(2x,i4)
           WRITE (50,1866) (theta(j),j=1,k),(sigf(jg), jg= 1,k),
     .  (kappa(jk), jk=1,K), (LAMBDA(JL),jl=1,K)
      WRITE (50,2066) (SIGMA(im,mj), mj = 1,M)
         do 2001 it = 0,T+1
 2001 write (50,2267) it,(z(it,nfa), nfa=1,k)
         go to 6466
 2066 format (6(1x,f13.10) )
 2267 format (i3,3(1x,f12.9) )
 1866 format (12(1x,f11.8) )
           end if
 7772   continue
 6466   continue
*****
 
        THETA(kneed) = DGCIN(prob,2,iptsig,table1,table2)
*       print 6686, kneed,theta(kneed)
 6686 format (i3,2x,'theta draw:',f12.8)
  1     CONTINUE
       RETURN
       END

       DOUBLE PRECISION FUNCTION PDFTH(x)
       IMPLICIT REAL * 8 (A-H,L,O-Z)
        DOUBLE PRECISION mnlpr,kvpri,kspri
        INTEGER T
       DOUBLE PRECISION KAPPA,klg,BSI(1000)
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /USE/ TWO,TOL
      COMMON /PRIORS/ tvpri(3),tspri(3),SIGPRV,
     c   SIGPRI(5,5),mnlpr(3),siglpr(3),vpris(3),spris(3),
     c   KVPRI(3),KSPRI(3)
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /SOME/ SUM
      COMMON /TCONS/ xtkons,zlnt,ntimth
      DIMENSION A(500,5,3),B(500,5,3)
      DIMENSION ep(500,5),facvar(3),t1(5)
         pi = 3.141592653589793238462643 d 0
       ntimth = ntimth + 1

      do 101 kf = 1,K
        FACVAR(kf) = sigf(kf) * sigf(kf)
        xkls = (kappa(kf) + LAMBDA(kf)) * (kappa(kf) + LAMBDA(kf))
        GAMMA(kf) = DSQRT(xkls + TWO * FACVAR(KF))
        klg = kappa(kf) + LAMBDA(KF) + GAMMA(KF)
        tg = two * gamma(kf)
       if (kf.eq.kneed) then
        aexp = (two * kappa(kf) * x ) / facvar(kf)
         else
        aexp = (two * kappa(kf) * theta(kf) ) / facvar(kf)
       end if
      DO 101 it = 1,T
      DO 101 ib = 1,M
       expo1 = gamma(kf)*smt(it,ib)
       expo2 = 0.5d 0 * klg * smt(it,ib)
        if (expo1.lt.600.0 d0.and.expo2.lt.600.0 d0) then
         ant = tg * dexp(expo2)
         adt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two * ( dexp(expo1) - 1.0 d0 )
         bdt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        else if (expo1.ge.600.0 d0.and.expo2.ge.600.0 d0) then
         ant = tg
         adt = klg * dexp(expo1 - expo2)
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two
         bdt = klg
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        else if (expo1.lt.600.0 d0.and.expo2.ge.600.0 d0) then
         ant = tg
         adt = klg *  dexp(expo1 - expo2)
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two * ( dexp(expo1) - 1.0 d0 )
         bdt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        else if (expo1.ge.600.0 d0.and.expo2.lt.600.0 d0) then
         ant = tg * dexp(expo2 - expo1)
         adt = klg
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two
         bdt = klg
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        end if
  101 CONTINUE

        s = 0.0 d 0

      do 3 it = 1,T
      DO 2, nb = 1,M
          pv = 0.0 d 0
      DO 1, nf = 1,K
         pv = pv + A(it,nb,nf) - B(it,nb,nf) * Z(it,nf)
  1   CONTINUE
         ep(it,nb) = r(it,nb) + pv
  2   CONTINUE
      do 22, j = 1,m
  22     t1(j) = 0.0 d 0
          do 31, ij = 1,M
          do 31, j = 1,M
 31       t1(ij) = t1(ij) + ep(it,j) *SIGINV(j,ij)
        do 32, j = 1,m
 32       s = s + t1(j) * ep(it,j)
  3     CONTINUE
        doug1 = -0.5 d0 * s
*  ft is the term from the likelihood in P(kappa(kneed) | .)
* posterior is likelihood * Prior
* Next evaluate p(z(k)|theta(k)) as fn of theta.

*       if (kneed.eq.2) then
*            smmm = .98 d 0
*       else 
*            smmm = 0.87 d 0
*       end if
        vf = sigf(kneed) * sigf(kneed)
        dt = 1.0 d 0 - DEXP(-kappa(kneed)*(1.0 d 0/52.14 d 0) )
        CIRC = (two * kappa(kneed)) / (vf * dt)
        CIRQ = ( (two*kappa(kneed)*x)/vf ) - 1.0 d 0
       ACIRQ = dabs(CIRQ)
 5555   pzdt = 1.0 d 0

       sdoug2 = 0.0 d 0
**  could put in the unconditional p(1,kneed) found in CIR ***
        do 3301 it = 2,T
      CIRU = CIRC * z(it-1,kneed) * DEXP(-kappa(kneed)*
     .    (1.0 d 0/52.14 d 0) )
      CIRV = CIRC * z(it,kneed)
      bbx = (CIRV / CIRU) ** (.5 d 0 * CIRQ)
      bx = two * ( (CIRU * CIRV)**.5 d0)
       if (bx.gt.1450.0 d 0) then
           zec = (8.0 d 0 * bx) * (8.0 d 0 * bx) * (8.0 d 0 * bx)
           BEST = 1.0 d 0 / dsqrt(two*pi*bx)
           bxmu = 4.0 d 0 * acirq*acirq
           best = best * (1.0 d 0 - (bxmu-1.0 d 0)/(8.0 d 0 * bx) +
     e ((bxmu-1.0 d 0)*(bxmu-9.0 d 0))/(two*(8.0 d 0*bx)*(8.0 d 0*bx))
     e - ((bxmu-1.0 d 0)*(bxmu-9.0 d 0)*(bxmu-25.0 d 0))/
     e  (6.0 d 0 *zec))
         if (best.lt.0.0 d 0) then
             best = 1.0 d -6
*            print 2398, best
          end if
 2398 format (' neg best approx used: ',f22.12)
      else
      if (ACIRQ.lt.1.0 d 0) then
         iu = 1
         go to 313
       end if
*  Note x must be positive, and xnu must fall between 0 and 1.  To allow
*  for nu values greater than 1, we have to take the integer amount of
*  draws.  The last one will be our boy.
       diu = dint(ACIRQ)
       iu = int(diu) + 1
       ACIRQ = DMOD(ACIRQ,diu)
 3136 format (' bx = ',f15.4)
 313   CALL DBSIES(ACIRQ,bx,iu,bsi)
* p. 10 in the special functions manual.
      BEST = BSI(iu)
        END IF
 4565 bby = CIRC * DEXP(-CIRU-CIRV+bx)
      st = (bby * bbx * BEST)
      doug2 = dlog(st)
      pzdt = pzdt * st * smmm
      sdoug2 = sdoug2 + doug2
 3301 CONTINUE
*       print 7389, x, sdoug2
 7389   format('at x: ',f10.5,' sdoug2: ',f30.10)
** now it is dexp(doug1+sdoug2)

      if (ntimth.eq.1) then 
          zlnt = (1.0 d 0 / pzdt) * .1 d -5
*         print 3456, pzdt
       end if
      if (ntimth.eq.1.and.kneed.eq.2) then
          zlnt = (1.0 d 0 / pzdt) * .1 d -5
*         print 3456, pzdt
 3456 format (' pzdt: ',f55.15)
       end if
 8760 format (' likelihood part: ',1(1x,f19.7) )
 8762 format (' pdf theta prts: ',3(1x,f19.7) )
      if (ntimth.eq.1) xtkons = doug1 + sdoug2
      PDFTH = dexp(doug1 + sdoug2 - xtkons)
      RETURN
      END
      SUBROUTINE DRWSIG
      IMPLICIT REAL * 8 (A-H,L,O-Z)
        DOUBLE PRECISION mnlpr,kappa,kvpri,kspri
       DOUBLE PRECISION LM(5,5),LMP(5,5),LMPQ(5,5)
      DIMENSION rc(1),U(5,5),varml(5,5)
     d,UP(5,5),V(5,5),VMI(5,5),UUP(5,5),eps(500,5),epspep(5,5)
        INTEGER T
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /PRIORS/ tvpri(3),tspri(3),SIGPRV,
     c   SIGPRI(5,5),mnlpr(3),siglpr(3),vpris(3),spris(3),
     c   KVPRI(3),KSPRI(3)
      COMMON /use/ two,tol
      COMMON /USFPAR/ A(500,5,3),B(500,5,3)
      DO 52 ir = 1,M
      DO 52 ic = 1,M
 52   varml(ir,ic) = 0.0 d 0
*     Form the T by M matrix of residuals.
      do 10 i = 1,t
      do 10 j = 1,M
        pv = 0.0 d 0
      do 1 jk = 1,K
  1     pv = pv + A(i,j,jk) - B(i,j,jk) * z(i,jk)
 10    eps(i,j) = r(i,j) + pv
      do 5 it = 1,t
      do 4 ir = 1,M
      do 4 ic = 1,M
      epspep(ir,ic) = eps(it,ir) * eps(it,ic)
  4   varml(ir,ic) = varml(ir,ic) + epspep(ir,ic)
  5   CONTINUE
      do 6 ir = 1,m
      do 6 ic = 1,m
  6   vmi(ir,ic) = varml(ir,ic)
      vt = dfloat(t)
        do 1099 jk = 1,M 
 1099   print 1096, (VMI(jk,jc), jc = 1,M)
 1096 format (5(1x,f26.19))
      CALL DLINRG(M,VMI,5,V,5)
      CALL DCHFAC(M,V,5,tol,irank,LM,5)
      do 7 i = 1,M
      do 7 j = 1,M
       LMP(i,j) = LM(j,i)
      if (i.lt.j) U(i,j) = 0.0 d 0
      if (i.gt.j) U(i,j) = drnnof()
      if (i.eq.j) then
          df = vt - dfloat(i)
          call drnchi(1,df,rc)
          U(i,j) = dsqrt(rc(1))
      end if
 7    CONTINUE
      do 12 i = 1,M
      do 12 j = 1,M
 12    UP(i,j) = U(j,i)
      do 13 i = 1,M
      do 13 j = 1,M
       uup(i,j) = 0.0 d 0
      do 13 kb = 1,M
  13  uup(i,j) = uup(i,j) + u(i,kb) * up(kb,j)
      do 14 i = 1,M
      do 14 j = 1,M
      LMPQ(i,j) = 0.0 d 0
      do 14 kb = 1,M
  14  LMPQ(i,j) = LMPQ(i,j) + LMP(I,kb) * uup(kb,j)
      do 15 i = 1,M
      do 15 j = 1,M
       siginv(i,j) = 0.0 d 0
      do 15 kb = 1,M
 15   SIGINV(i,j) = SIGINV(i,j) + LMPQ(i,kb) * LM(kb,j)
*       do 1299 jk = 1,M 
*1299   print 1096, (SIGINV(jk,jc), jc = 1,M)
      CALL DLINRG(M,SIGINV,5,SIGMA,5)
      RETURN
      END
        SUBROUTINE GETK
        IMPLICIT REAL * 8 (A-H,L,O-Z)
        DIMENSION TABLE1(500),TABLE2(500),xv(0:10000),pval(10000),
     d   udifg(900)
        DOUBLE PRECISION mnlpr,kappa,kvpri,kspri
        INTEGER T  
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /PRIORS/ tvpri(3),tspri(3),SIGPRV,
     c   SIGPRI(5,5),mnlpr(3),siglpr(3),vpris(3),spris(3),
     c   KVPRI(3),KSPRI(3)
      COMMON /USE/ TWO,TOL
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /KCONS/ xkkons,ntimth
       EXTERNAL pdfk
       diffi = 0.07 d -2
       iptsig = 90
        errabs = 0.0 d 0
        zero = 0.0 d 0
        errel = .1 d -1
       do 1 kneed = 1,K
          diffo = kappa(kneed) * .05 d -1
*         if (diffo.lt.diffi) diffo = diffi
          if (kneed.eq.3) diffo = .001
          pvmax = -1.0 d -6
 2347      ntimth = 0
        garb = pdfk(kappa(kneed))
          xv(0) = 0.92 d 0 * kappa(kneed) 
*         if (kneed.eq.1) xv(0) = .8 * kappa(kneed)
*         JNEED = 2.2 d 0 * (kappa(kneed) - xv(0) ) / diffo
          JNEED = 10000
 4565      if (xv(0).lt.zero) xv(0) = zero
          pvmax = -1.0 d -6
         DO 1010 jt = 1,JNEED
          xv(jt) = xv(jt-1) + diffo
         pval(jt) = pdfk(xv(jt))
*         print 6758, xv(jt),pval(jt)
 6758 format (1x,f9.6,1x,f29.2)
            tpv = pvmax - .1 d -4
            if (pval(jt).gt.tpv) then 
               oldmax = pvmax
               pvmax = pval(jt)
               jmax = jt
               xmxa = pvmax * 0.0005 d 0
 765  format (2(1x,f52.43) )
             else
               if (pval(jt).lt.xmxa) go to 1011
             END IF
 1010  CONTINUE
 1011  trylas = xv(jt)
*      print 6378
 6378  format(' I got to this point ')
       table1(iptsig) = trylas
       table2(iptsig) = pval(jt)
         DO 102 iu = 1,jt
            if (pval(iu).gt.xmxa) then
               if (iu.eq.1) then
               print 7865
 7865 format (' starting point was too large.')
                  xv(0) = .8 d 0 * xv(0)
                  GO TO 4565
               end if
              go to 1021
             end if
  102    CONTINUE
       print 6378
 1021  trystr = xv(iu-1)
       table2(1) = pval(iu-1)
       iptsig = 90
       iepts = 78
       difg = (trylas - trystr) / dfloat(iptsig + iepts)
        if (difg.le.zero) then
           difg = dabs(difg)
*           print 1026, trylas,trystr,difg
        end if
 1026 format (' negative difg. trylas: ',f12.8,' trystr: ',f12.7,' difg:
     f   ',f9.6)
 8000  do 8907 itp = 1,iptsig
 8907 udifg(itp) = difg
       table1(1) = trystr
          udifg(2) = 20.0 d 0 * difg
          udifg(3) = 10.0 d 0 * difg
          udifg(4) = 6.0 d 0 * difg
          udifg(5) = 3.0 d 0 * difg
          udifg(iptsig-1) = udifg(2)
          udifg(iptsig-2) = udifg(3)
          udifg(iptsig-3) = udifg(4)
          udifg(iptsig-4) = udifg(5)
       do 76 ig = 2,iptsig
       table1(ig) = table1(ig-1) + udifg(ig)
       table2(ig) = pdfk(table1(ig))
*       PRINT 7676, TABLE1(IG),TABLE2(IG) 
 7676 format (2(1x,f22.6) )
  76   CONTINUE
       prob = drnunf()

***** 
        do 7772, iv = 1, iptsig - 1
         if(table1(iv+1).le.table1(iv)) then
           write(50, 7990)
 7990   format(' problem in kappa ')
           write(50, 7991) kneed
 7991   format(2x,i4)
           WRITE (50,1866) (theta(j),j=1,k),(sigf(jg), jg= 1,k),
     .  (kappa(jk), jk=1,K), (LAMBDA(JL),jl=1,K)
      WRITE (50,2066) (SIGMA(im,mj), mj = 1,M)
         do 2001 it = 0,T+1
 2001 write (50,2267) it,(z(it,nfa), nfa=1,k)
         go to 6466
 2066 format (6(1x,f13.10) )
 2267 format (i3,3(1x,f12.9) )
 1866 format (12(1x,f11.8) )
           end if
 7772   continue
 6466   continue
*****

        KAPPA(kneed) = DGCIN(prob,2,iptsig,table1,table2)
  1     CONTINUE
        RETURN
        END
      DOUBLE PRECISION FUNCTION PDFK(x)
       IMPLICIT REAL * 8 (A-H,L,O-Z)
        DOUBLE PRECISION mnlpr,kvpri,kspri
        INTEGER T
       DOUBLE PRECISION KAPPA,klg,bsi(1000)
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /USE/ TWO,TOL
      COMMON /PRIORS/ tvpri(3),tspri(3),SIGPRV,
     c   SIGPRI(5,5),mnlpr(3),siglpr(3),vpris(3),spris(3),
     c   KVPRI(3),KSPRI(3)
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /SOME/ SUM
      COMMON /KCONS/ xkkons,ntimth
      DIMENSION A(500,5,3),B(500,5,3)
      DIMENSION ep(500,5),facvar(3),t1(5)
         pi = 3.141592653589793238462643 d 0
*       print 6545, x
 6545 format (' inside pdfk x = ',f12.7)
       ntimth = ntimth + 1
      do 101 kf = 1,K
        FACVAR(kf) = sigf(kf) * sigf(kf)
       if (kf.eq.kneed) then
         xkls = (x + LAMBDA(kf)) * (x + LAMBDA(kf))
         GAMMA(kf) = DSQRT(xkls + TWO * FACVAR(KF))
        klg = x + LAMBDA(KF) + GAMMA(KF)
        aexp = (two * x * theta(kf) ) / facvar(kf)
         else
          xkls = (kappa(kf) + LAMBDA(kf)) * (kappa(kf) + LAMBDA(kf) )
        GAMMA(kf) = DSQRT(xkls + TWO * FACVAR(KF))
        klg = kappa(kf) + LAMBDA(KF) + GAMMA(KF)
        aexp = (two * kappa(kf) * theta(kf) ) / facvar(kf)
         end if
        tg = two * gamma(kf)
      DO 101 it = 1,T
      DO 101 ib = 1,M
       expo1 = gamma(kf)*smt(it,ib)
       expo2 = 0.5d 0 * klg * smt(it,ib)
        if (expo1.lt.600.0 d0.and.expo2.lt.600.0 d0) then
         ant = tg * dexp(expo2) 
         adt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two * ( dexp(expo1) - 1.0 d0 )
         bdt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        else if (expo1.ge.600.0 d0.and.expo2.ge.600.0 d0) then
         ant = tg
         adt = klg * dexp(expo1 - expo2)
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two
         bdt = klg
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        else if (expo1.lt.600.0 d0.and.expo2.ge.600.0 d0) then
         ant = tg
         adt = klg *  dexp(expo1 - expo2)
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two * ( dexp(expo1) - 1.0 d0 )
         bdt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        else if (expo1.ge.600.0 d0.and.expo2.lt.600.0 d0) then
         ant = tg * dexp(expo2 - expo1)
         adt = klg
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two
         bdt = klg
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        end if
  101 CONTINUE
         s = 0.0 d 0
       do 3 it = 1,T
      DO 2, nb = 1,M
          pv = 0.0 d 0
      DO 1, nf = 1,K
         pv = pv + A(it,nb,nf) - B(it,nb,nf) * Z(it,nf)
  1   CONTINUE
         ep(it,nb) = r(it,nb) + pv
  2   CONTINUE
      do 22, j = 1,m 
  22     t1(j) = 0.0 d 0
          do 31, ij = 1,M
          do 31, j = 1,M
 31       t1(ij) = t1(ij) + ep(it,j) *SIGINV(j,ij)  
        do 32, j = 1,m
 32       s = s + t1(j) * ep(it,j)
  3     CONTINUE
*      if (ntimth.eq.1) xkkons = .5 d 0 * s
*      print 4564, ntimth,xkkons
 4564 format (' ntimth: ',i4,' xkkons: ',f40.11)
*     ft = dexp(-.5 * s + xkkons)
       doug1 = -.5 d 0 * s
*      WRITE (16,36) x,ft
 36    format (' pdfk x: ',f12.7,' ft = ',f12.6)
*  ft is the term from the likelihood in P(kappa(kneed) | .)
* posterior is likelihood * Prior
* Next evaluate p(z(k)|kappa(k)) as fn of kappa.
* 
        pzdk = 1.0 d 0
        vf = sigf(kneed) * sigf(kneed)
        dt = 1.0 d 0 - DEXP(-x*(1.0 d 0/52.14 d 0) )
        CIRC = (two * x) / (vf * dt)
        CIRQ = ( (two*x*theta(kneed))/vf ) - 1.0 d 0
       ACIRQ = dabs(CIRQ)
        sdoug2 = 0.0 d 0
        do 3301 it = 2,T
        CIRU = CIRC * z(it-1,kneed) * DEXP(-x*(1.0 d 0/52.14 d 0) )
        CIRV = CIRC * z(it,kneed)
      bbx = (CIRV / CIRU) ** (.5 d 0 * CIRQ)
      bx = two * ( (CIRU * CIRV)**.5 d0)
       if (bx.gt.1450.0 d 0) then
           zec = (8.0 d 0 * bx) * (8.0 d 0 * bx) * (8.0 d 0 * bx)
           BEST = 1.0 d 0 / dsqrt(two*pi*bx)
           bxmu = 4.0 d 0 * acirq*acirq
           best = best * (1.0 d 0 - (bxmu-1.0 d 0)/(8.0 d 0 * bx) +
     e ((bxmu-1.0 d 0)*(bxmu-9.0 d 0))/(two*(8.0 d 0*bx)*(8.0 d 0*bx))
     e - ((bxmu-1.0 d 0)*(bxmu-9.0 d 0)*(bxmu-25.0 d 0))/
     e  (6.0 d 0 *zec))
*     print 6413, best
 6413 format (' from approx: best = ',f22.5)
         if (best.lt.0.0 d 0) then
             best = 1.0 d -6
*            print 2398, best
          end if
 2398 format (' neg best approx used: ',f22.12)
      else
      if (ACIRQ.lt.1.0 d 0) then
         iu = 1
         go to 313
       end if
*     do 310 i = 2,200
*        if (dfloat(i).gt.ACIRQ.and.dfloat(i-1).lt.ACIRQ ) go to 312
*310   CONTINUE
*     Print 3166, ACIRQ
 3166 format (' problem - ACIRQ = ',f12.5)
*  Note x must be positive, and xnu must fall between 0 and 1.  To allow
*  for nu values greater than 1, we have to take the integer amount of
*  draws.  The last one will be our boy.
* 312  iu = i-1
*      ACIRQ = ACIRQ - dfloat(iu)
       diu = dint(ACIRQ)
       iu = int(diu) + 1
       ACIRQ = DMOD(ACIRQ,diu)
 3136 format (' bx = ',f15.4)
 313   CALL DBSIES(ACIRQ,bx,iu,bsi)
* p. 10 in the special functions manual.
      BEST = BSI(iu)
*        print 6543, best
 6543 format (' from IMSL: best = ',f22.5)
       END IF
*      print 7393, it,CIRC,CIRU,CIRV,bx
 7393  format(' it = ',i2,1x,f30.5,1x,3(f16.8,1x))
 4565 bby = CIRC * DEXP(-CIRU-CIRV+bx)
      st = (bby * bbx * BEST)
       doug2 = dlog(st)
       sdoug2 = sdoug2 + doug2
*       smmm = 1.0 d 0
*       if (kneed.eq.2) smmm = 0.2 d -2
*     pzdk = pzdk * st * smmm
 3301 CONTINUE
*     if (ntimth.eq.1) then
*           zlnk = (1.0 d 0 / pzdk)
*         print 3456, pzdk
*       end if
 3456 format (' pzdk: ',f55.15)
*       pzdk = pzdk * zlnk
        if (ntimth.eq.1) xkkons = doug1 + sdoug2
       PDFK = dexp(doug1 + sdoug2 - xkkons)
*       print 654, x,pdfk
 654   format (f12.7,1x,f40.5)
      RETURN
      END
        SUBROUTINE GETL
        IMPLICIT REAL * 8 (A-H,L,O-Z)
         DIMENSION TABLE1(500),TABLE2(500),UDIFG(999)
        DOUBLE PRECISION mnlpr,kappa,kvpri,kspri,xv(0:10000),pval(10000)
        INTEGER T  
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /PRIORS/ tvpri(3),tspri(3),SIGPRV,
     c   SIGPRI(5,5),mnlpr(3),siglpr(3),vpris(3),spris(3),
     c   KVPRI(3),KSPRI(3)
      COMMON /USE/ TWO,TOL
      COMMON /CNST/ XMCONS,nocall
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /SOME/ SUM
      EXTERNAL DRAWL,pdfl
       errabs = 0.0 d 0
       zero = 0.0 d 0
       errel = .01 d 0
       diffo = .75 d -3
       iptsig = 90
       iepts = 78
c  kneed stands for the factor whose parameter we need to draw
       do 1, kneed = 1,K
        nocall = 0
        garb = pdfl(lambda(kneed))
        XMCONS = XMCONS * .73 d 0
        ndr = 0
          pvmax = -1.0 d -8
           xmxa = -99999.9 d 0
          xv(0) = lambda(kneed) - .027 d 0
         if (xv(0).gt.zero) xv(0) = -.01 d 0
*        JNEED = 0.09 d 0/ diffo
         JNEED = 10000
         DO 1010 jt = 1,JNEED
          xv(jt) = xv(jt-1) + diffo
         pval(jt) = pdfl(xv(jt))
*      print 1088, xv(jt),pval(jt)
 1088 format (f12.9,1x,f29.1)
            tpv = pvmax - .1 d -9
            if (pval(jt).gt.tpv) then
               oldmax = pvmax
               pvmax = pval(jt)
               jmax = jt
               xmxa = pvmax * 0.0005 d 0
 765  format (2(1x,f52.43) )
             else
*            print 6758, xmxa
 6758 format (' xmxa: ',f17.9)
               if (pval(jt).lt.xmxa) go to 1011
             end if
 1010  CONTINUE
 1011  trylas = xv(jt)
       table1(iptsig) = trylas
       table2(iptsig) = pval(jt)
         DO 102 iu = 1,jt
            if (pval(iu).gt.xmxa) go to 1021
  102    CONTINUE
 1021  trystr = xv(iu-1)
       table2(1) = pval(iu-1)
       difg = (trylas - trystr) / dfloat(iptsig + iepts)
 8000  do 8907 itp = 1,iptsig
 8907 udifg(itp) = difg
       table1(1) = trystr
          udifg(2) = 20.0 d 0 * difg
          udifg(3) = 10.0 d 0 * difg
          udifg(4) = 6.0 d 0 * difg
          udifg(5) = 3.0 d 0 * difg
          udifg(iptsig-1) = udifg(2)
          udifg(iptsig-2) = udifg(3)
          udifg(iptsig-3) = udifg(4)
          udifg(iptsig-4) = udifg(5)
       do 76 ig = 2,iptsig
       table1(ig) = table1(ig-1) + udifg(ig)
  76   table2(ig) = pdfl(table1(ig))
*      do 7622 iup = 1,iptsig
*7622 print 8765, iup,table1(iup),table2(iup)
 8765 format (1x,i3,1x,f10.7,1x,f70.12 )
       prob = drnunf()

***** 
        do 7772, iv = 1, iptsig - 1
         if(table1(iv+1).le.table1(iv)) then
           write(50, 7990)
 7990   format(' problem in lambda ')
           write(50, 7991) kneed
 7991   format(2x,i4)
           WRITE (50,1866) (theta(j),j=1,k),(sigf(jg), jg= 1,k),
     .  (kappa(jk), jk=1,K), (LAMBDA(JL),jl=1,K)
      WRITE (50,2066) (SIGMA(im,mj), mj = 1,M)
         do 2001 it = 0,T+1
 2001 write (50,2267) it,(z(it,nfa), nfa=1,k)
         go to 6466
 2066 format (6(1x,f13.10) )
 2267 format (i3,3(1x,f12.9) )
 1866 format (12(1x,f11.8) )
           end if
 7772   continue
 6466   continue
*****

        LAMBDA(kneed) = DGCIN(prob,2,iptsig,table1,table2)
 1      CONTINUE
        RETURN
        END
      DOUBLE PRECISION FUNCTION PDFL(x)
       IMPLICIT REAL * 8 (A-H,L,O-Z)
       INTEGER T
        DIMENSION T1(5),facvar(3),ep(500,5)
        DOUBLE PRECISION mnlpr,kappa,kvpri,kspri
       DOUBLE PRECISION lambda,klg
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /USE/ TWO,TOL
      DIMENSION A(500,5,3),B(500,5,3)
      COMMON /PRIORS/ tvpri(3),tspri(3),SIGPRV,
     c   SIGPRI(5,5),mnlpr(3),siglpr(3),vpris(3),spris(3),
     c   KVPRI(3),KSPRI(3)
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /SOME/ SUM
      COMMON /CNST/ XMCONS,nocall
       nocall = nocall + 1
      do 101 kf = 1,K
        FACVAR(kf) = sigf(kf) * sigf(kf)
       if (kf.eq.kneed) then
            gt = (KAPPA(kf) + x) * (KAPPA(kf) + x)
         GAMMA(kf) = DSQRT(gt + TWO * FACVAR(KF))
        klg = kAPPA(kf) + x + GAMMA(KF)
        aexp = (two * KAPPA(kf) * theta(kf) ) / facvar(kf)
         else
         gt = (KAPPA(kf) + LAMBDA(kf) ) * (KAPPA(kf) + LAMBDA(kf) )
        GAMMA(kf) = DSQRT(gt + TWO * FACVAR(KF))
        klg = kappa(kf) + LAMBDA(KF) + GAMMA(KF)
        aexp = (two * kappa(kf) * theta(kf) ) / facvar(kf)
         end if
        tg = two * gamma(kf)
      DO 101 it = 1,T
      DO 101 ib = 1,M
       expo1 = gamma(kf)*smt(it,ib)
       expo2 = 0.5d 0 * klg * smt(it,ib)
        if (expo1.lt.600.0 d0.and.expo2.lt.600.0 d0) then
         ant = tg * dexp(expo2)
         adt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two * ( dexp(expo1) - 1.0 d0 )
         bdt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        else if (expo1.ge.600.0 d0.and.expo2.ge.600.0 d0) then
         ant = tg
         adt = klg * dexp(expo1 - expo2)
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two
         bdt = klg
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        else if (expo1.lt.600.0 d0.and.expo2.ge.600.0 d0) then
         ant = tg
         adt = klg *  dexp(expo1 - expo2)
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two * ( dexp(expo1) - 1.0 d0 )
         bdt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        else if (expo1.ge.600.0 d0.and.expo2.lt.600.0 d0) then
         ant = tg * dexp(expo2 - expo1)
         adt = klg
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two
         bdt = klg
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
        end if
  101 CONTINUE
         s = 0.0 d 0
       do 3 it = 1,T
      DO 2, nb = 1,M
       pv = 0.0 d 0
      DO 1, nf = 1,K
         pv = pv + A(it,nb,nf) - B(it,nb,nf) * Z(it,nf)
  1   CONTINUE
        ep(it,nb) = R(it,nb) + pv
  2   CONTINUE
*        write (16,1654) it,(ep(it,jb), jb=1,M)
 1654 format (1x,i4,4(1x,f22.13) )
      do 22, j = 1,m
  22     t1(j) = 0.0 d 0
          do 31, ij = 1,M
          do 31, j = 1,M
 31       t1(ij) = t1(ij) + ep(it,j) *SIGINV(j,ij)
        do 32, j = 1,m
 32       s = s + t1(j) * ep(it,j)
*      WRITE (16,2654) it,s
 2654 format (' S: ',i3,1x,f26.11)
  3     CONTINUE
      if (nocall.eq.1) xmcons = s * 0.81 d 0
      ft = dexp(-.5 * s + xmcons)
*  ft is the term from the likelihood in P(lambda(kneed) | .)
      pdfl = ft
*         if (pdfl.gt.0.0 d 0.or.pdfl.lt.9999.9 d 0) go to 987
*         PRINT 765, x,s,st,pdfl
 765   format (4(1x,f40.30) )
 987   RETURN
       END
        SUBROUTINE GETSF
        IMPLICIT REAL * 8 (A-H,L,O-Z)
        DIMENSION TABLE1(500),TABLE2(500)
        DOUBLE PRECISION mnlpr,kappa,kspri,kvpri
        DOUBLE PRECISION xv(0:50000),pval(50000),UDIFG(999)
        INTEGER T  
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /PRIORS/ tvpri(3),tspri(3),SIGPRV,
     c   SIGPRI(5,5),mnlpr(3),siglpr(3),vpris(3),spris(3),
     c   KVPRI(3),KSPRI(3)
      COMMON /USE/ TWO,TOL
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /SCONS/ xsskons,ncl
      COMMON /SOME/ SUM
      EXTERNAL DRAWS,pdfs
        errabs = 0.0 d 0
        zero = 0.0 d 0
        errel = .01 d 0
        addt = 0.032 d 0
       do 1 kneed = 1,K
         ncl = 0
         gar = pdfs(sigf(kneed) )
         if (kneed.eq.2) then
             diffo = .023 d -2
             addt = 0.033 d 0
         end if
         if (kneed.eq.3) then
               diffo = .5 d -2
               addt = 0.012
           end if
         if (kneed.eq.1) then
             diffo = 0.6 d -2
             addt = .36 d 0 * sigf(kneed)
          end if
        icntfn = 0
*   Set up the table
 2347     pvmax = -10.0 d 0
           xmxa = -99999.9 d 0
*          JNEED = (2.0 d 0 * addt) / diffo
           JNEED = 10000
          xv(0) = sigf(kneed) - addt
          xv(0) = dmax1(xv(0),0.0 d 0)
         DO 1010 jt = 1,JNEED + 1000
          xv(jt) = xv(jt-1) + diffo
          if (xv(jt).gt.1.2 d 0) then
                 addt = addt * 1.7 d 0
                 diffo = diffo * 0.5
                 go to 2347
             end if
         pval(jt) = pdfs(xv(jt))
*        print 765, xv(jt),pval(jt)
            tpv = pvmax - .1 d -23
            if (pval(jt).gt.tpv) then
               oldmax = pvmax
               pvmax = pval(jt)
               jmax = jt
               xmxa = pvmax * 0.0005 d 0
 765  format (1x,f8.6,1x,f45.33)
             else
               if (pval(jt).lt.xmxa) go to 1011
             end if
        if (jt.gt.49990) then
           icntfn = icntfn + 1
           print 8979
 8979 format (' seg fault warning in sigf.')
        addt = 0.025 / dfloat(icntfn)
        diffo = .2 d -5 / dfloat(icntfn)
         if (icntfn.gt.6) stop 7654
        go to 2347
        end if
 1010  CONTINUE
 1011  trylas = xv(jt)
         DO 102 iu = 1,jt
            if (pval(iu).gt.xmxa) go to 1021
  102    CONTINUE
 1021  trystr = xv(iu-1)
       table2(1) = pval(iu-1)
       if (trylas.le.trystr) then
          diffo = diffo / 1.08 d 0
          addt = 1.03 d 0 * addt
            BIGNO = BIGNO + 25.0 d 0
          go to 2347
        end if
       iptsig = 60
       iepts = 78
       difg = (trylas - trystr) / dfloat(iptsig + iepts)
 8000  do 8907 itp = 1,iptsig
 8907 udifg(itp) = difg
       table1(1) = trystr
          udifg(2) = 20.0 d 0 * difg
          udifg(3) = 10.0 d 0 * difg
          udifg(4) = 6.0 d 0 * difg
          udifg(5) = 3.0 d 0 * difg
          udifg(iptsig-1) = udifg(2)
          udifg(iptsig-2) = udifg(3)
          udifg(iptsig-3) = udifg(4)
          udifg(iptsig-4) = udifg(5)
       do 76 ig = 2,iptsig
       table1(ig) = table1(ig-1) + udifg(ig)
  76   table2(ig) = pdfs(table1(ig))
       prob = drnunf()
*        do 276 io= 1,iptsig
*276     print 2766, io,table1(io),table2(io)
 2766 format (i4,21x,f8.6,1x,f66.3)
*        print 678, prob
 678     format (' unif draw: ',f9.7)

***** 
        do 7772, iv = 1, iptsig - 1
         if(table1(iv+1).le.table1(iv)) then
           write(50, 7990)
 7990   format(' problem in sigf ')
           write(50, 7991) kneed
 7991   format(2x,i4)
           WRITE (50,1866) (theta(j),j=1,k),(sigf(jg), jg= 1,k),
     .  (kappa(jk), jk=1,K), (LAMBDA(JL),jl=1,K)
      WRITE (50,2066) (SIGMA(im,mj), mj = 1,M)
         do 2001 it = 0,T+1
 2001 write (50,2267) it,(z(it,nfa), nfa=1,k)
         go to 6466
 2066 format (6(1x,f13.10) )
 2267 format (i3,3(1x,f12.9) )
 1866 format (12(1x,f11.8) )
           end if
 7772   continue
 6466   continue
*****

        SIGF(kneed) = DGCIN(prob,2,iptsig,table1,table2)
 1      CONTINUE
        RETURN
        END
      DOUBLE PRECISION FUNCTION PDFS(x)
       IMPLICIT REAL * 8 (A-H,L,O-Z)
        INTEGER T
        DIMENSION t1(6),ep(400,6),facvar(3)
        DOUBLE PRECISION mnlpr,kvpri,kspri,bsi(1000)
       DOUBLE PRECISION KAPPA,klg
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /PRIORS/ tspri(3),tvpri(3),SIGPRV,
     c   SIGPRI(5,5),mnlpr(3),siglpr(3),vpris(3),spris(3),
     c   KVPRI(3),KSPRI(3)
      COMMON /USE/ TWO,TOL
      DIMENSION A(500,5,3),B(500,5,3)
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /SCONS/ xsskons,ncl
      COMMON /SOME/ SUM
        pi = 3.141592653589793238462643 d 0
*       print 2796, theta(1),theta(2)
*       print 2796, kappa(1),kappa(2)
*       print 2796, lambda(1),lambda(2)
*       print 2796, sigf(1),sigf(2)
 2796  format(2(2x,f20.10))
        ncl = ncl + 1
      do 101 kf = 1,K
        if (kf.eq.kneed) then
           facvar(kf) = x * x
         ELSE
        FACVAR(kf) = sigf(kf) * sigf(kf)
         end if
          gt = (KAPPA(kf) + LAMBDA(kf) ) * (KAPPA(kf) + LAMBDA(kf) )
         GAMMA(kf) = DSQRT(GT + TWO * FACVAR(KF))
        klg = kAPPA(kf) + LAMBDA(kf) + GAMMA(KF)
        aexp = (two * KAPPA(kf) * theta(kf) ) / facvar(kf)
        tg = two * gamma(kf)
      DO 101 it = 1,T
      DO 101 ib = 1,M
       expo1 = gamma(kf)*smt(it,ib)
       expo2 = 0.5d 0 * klg * smt(it,ib)
         ant = tg * dexp(expo2)
         adt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         A(it,ib,kf) = ( aexp * dlog(ant/adt) ) / smt(it,ib)
         bnt = two * ( dexp(expo1) - 1.0 d0 )
         bdt = tg + klg * ( dexp(expo1) - 1.0 d0 )
         B(it,ib,kf) = ( bnt/bdt ) / smt(it,ib)
  101 CONTINUE
         s = 0.0 d 0
       do 3 it = 1,T
      DO 2, nb = 1,M
         pv = 0.0 d 0
      DO 1, nf = 1,K
         pv = pv + A(it,nb,nf) - B(it,nb,nf) * Z(it,nf)
  1   CONTINUE
         ep(it,nb) = R(it,nb) + pv
  2   CONTINUE
      do 22, j = 1,m
  22     t1(j) = 0.0 d 0
          do 31, ij = 1,M
          do 31, j = 1,M
 31       t1(ij) = t1(ij) + ep(it,j) *SIGINV(j,ij)
        do 32, j = 1,m
 32       s = s + t1(j) * ep(it,j)
  3     CONTINUE
*       IF (NCL.eq.1) then
*            BIGNO = -1.00 d 0 * (1.0 d 0 - (0.5 d 0 * s))
*            print 3336, bigno
*3336 format (' bigno = ',f33.9)
*       END IF
*     ft = dexp(-.5 * s + BIGNO)

**new 7 13 98
      doug1 = -0.5 d0 * s
** end new

*  ft is the term from the likelihood in P(sigma(kneed) | .)
* posterior is likelihood * Prior
        pzds = 1.0 d 0
        vf = x * x
        dt = 1.0 d 0 - DEXP(-kappa(kneed) *(1.0 d 0/52.14 d 0) )
        CIRC = (two * kappa(kneed) ) / (vf * dt)
*       print 5412, two, kappa(kneed), vf, dt
 5412   format(4(f10.5))
        CIRQ = ( (two*kappa(kneed)*theta(kneed))/vf ) - 1.0 d 0
       ACIRQ = dabs(CIRQ)
        sdoug2 = 0.0 d 0
        do 3301 it = 2,T
        CIRU = CIRC * z(it-1,kneed) * DEXP(-kappa(kneed)
     p                 * (1.0 d 0/52.14 d 0) )
        CIRV = CIRC * z(it,kneed)
      bbx = (CIRV / CIRU) ** (.5 d 0 * CIRQ)
      bx = two * ( (CIRU * CIRV)**.5 d0)
       if (bx.gt.1450.0 d 0) then
           zec = (8.0 d 0 * bx) * (8.0 d 0 * bx) * (8.0 d 0 * bx)
           BEST = 1.0 d 0 / dsqrt(two*pi*bx)
           bxmu = 4.0 d 0 * acirq*acirq
           best = best * (1.0 d 0 - (bxmu-1.0 d 0)/(8.0 d 0 * bx) +
     e ((bxmu-1.0 d 0)*(bxmu-9.0 d 0))/(two*(8.0 d 0*bx)*(8.0 d 0*bx))
     e - ((bxmu-1.0 d 0)*(bxmu-9.0 d 0)*(bxmu-25.0 d 0))/
     e  (6.0 d 0 *zec))
*     print 6413, best
 6413 format (' from approx: best = ',f22.5)
         if (best.lt.0.0 d 0) then
             best = 1.0 d -6
*            print 2398, best
          end if
 2398 format (' neg best approx used: ',f22.12)
      else
      if (ACIRQ.lt.1.0 d 0) then
         iu = 1
         go to 313
       end if
*     do 310 i = 2,200
*        if (dfloat(i).gt.ACIRQ.and.dfloat(i-1).lt.ACIRQ ) go to 312
*310   CONTINUE
*     Print 3166, ACIRQ
 3166 format (' problem - ACIRQ = ',f12.5)
*  Note x must be positive, and xnu must fall between 0 and 1.  To allow
*  for nu values greater than 1, we have to take the integer amount of
*  draws.  The last one will be our boy.
* 312  iu = i-1
*      ACIRQ = ACIRQ - dfloat(iu)
       diu = dint(ACIRQ)
       iu = int(diu) + 1
       ACIRQ = DMOD(ACIRQ,diu)
 3136 format (' bx = ',f15.4)
 313   CALL DBSIES(ACIRQ,bx,iu,bsi)
* p. 10 in the special functions manual.
       BEST = BSI(iu)
*        print 6543, best
 6543 format (' from IMSL: best = ',f22.5)
       END IF
*      print 7393, it,CIRC,CIRU,CIRV,bx
 7393  format(' it = ',i2,1x,f30.5,1x,3(f16.8,1x))
       dbby1 = dlog(CIRC) + (-CIRU-CIRV+bx)
       dbbx1 = dlog(bbx)
       dbbz1 = dlog(BEST)
       doug2 = dbby1 + dbbx1 + dbbz1
       sdoug2 = sdoug2 + doug2
*      print 7399, it,dbby1,dbbx1,dbbz1
 7399  format(' it = ',i2,1x,3(f16.8,1x))
*       smmm = 1.0 d 0
*       if (kneed.eq.2) smmm = 0.2 d -2
*     pzdk = pzdk * st * smmm
 3301 CONTINUE
        if (ncl.eq.1) xsskons = (doug1 + sdoug2) * 1.0 d 0
        if (kneed.eq.2.and.ncl.eq.1) xsskons =
     e doug1 + sdoug2
        PDFS = dexp(doug1 + sdoug2 - xsskons)
*      print 7363, doug1, sdoug2,xsskons
 7363  format(4(2x,f15.10))
** end new
       RETURN
       END

        SUBROUTINE AUGZ
        IMPLICIT REAL * 8 (A-H,L,O-Z)
        DOUBLE PRECISION kappa
         DIMENSION TABLE1(500),TABLE2(500),xv(0:10000),pval(10000),
     d     udifg(999)
        INTEGER T
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /USE/ TWO,TOL
      COMMON /USFPAR/ A(500,5,3),B(500,5,3)
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /SOME/ SUM
      COMMON /LOWBOU/ zlow
      EXTERNAL DRAWZ0,pdfz0,DRAWZP,PDFZP
         diffo = 0.3 d -3
        errabs = 0.0 d 0
        zero = 0.0 d 0
        errel = .01 d 0
       do 1 kneed = 1,K
       azlb = z(0,kneed) - .4 d 0 * z(0,kneed)
*   Set up the table
2347      pvmax = -0.1 d -9
           xmxa = -99999.9 d 0
*         JNEED = .18 d 0 / diffo
          JNEED = 10000
          xv(0) = azlb
         DO 1210 jt = 1,JNEED
          xv(jt) = xv(jt-1) + diffo
         pval(jt) = pdfz0(xv(jt))
*       print 2366, xv(jt),pval(jt)
 2366 format (2f18.9)
            tpv = pvmax - .1 d -8
            if (pval(jt).gt.tpv) then
               oldmax = pvmax
               pvmax = pval(jt)
               xmxa = pvmax * 0.0005 d 0
 765  format (2(1x,f52.43) )
             else
               if (pval(jt).lt.xmxa) go to 1211
             end if
 1210  CONTINUE
 1211  trylas = xv(jt)
         DO 1022 iu = 1,jt
            if (pval(iu).gt.xmxa) go to 1221
 1022    CONTINUE
 1221  trystr = xv(iu-1)
*      print 6765, trystr,pval(iu),pvmax,trylas,pval(jt)
 6765 format (5f17.9)
       iptsig = 50
       iepts = 78
       difg = (trylas - trystr) / dfloat(iptsig + iepts)
       do 8907 itp = 1,iptsig
 8907 udifg(itp) = difg
       table1(1) = trystr
       table2(1) = pval(iu-1)
          udifg(2) = 20.0 d 0 * difg
          udifg(3) = 10.0 d 0 * difg
          udifg(4) = 6.0 d 0 * difg
          udifg(5) = 3.0 d 0 * difg
          udifg(49) = udifg(2)
          udifg(48) = udifg(3)
          udifg(47) = udifg(4)
          udifg(46) = udifg(5)
       do 76 ig = 2,iptsig
       table1(ig) = table1(ig-1) + udifg(ig)
  76   table2(ig) = pdfz0(table1(ig))
       prob = drnunf()

***** 
        do 7772, iv = 1, iptsig - 1
         if(table1(iv+1).le.table1(iv)) then
           write(50, 7990)
 7990   format(' problem in zee zero ')
           write(50, 7991) kneed
 7991   format(2x,i4)
           WRITE (50,1866) (theta(j),j=1,k),(sigf(jg), jg= 1,k),
     .  (kappa(jk), jk=1,K), (LAMBDA(JL),jl=1,K)
      WRITE (50,2066) (SIGMA(im,mj), mj = 1,M)
         do 2001 it = 0,T+1
 2001 write (50,2267) it,(z(it,nfa), nfa=1,k)
         go to 6466
 2066 format (6(1x,f13.10) )
 2267 format (i3,3(1x,f12.9) )
 1866 format (12(1x,f11.8) )
           end if
 7772   continue
 6466   continue
*****

        z(0,kneed) = DGCIN(prob,2,iptsig,table1,table2)
* Now get the T + 1 Z.
        diffo = 0.1 d -3
        azlb = z(T+1,kneed) - .2 d 0 * z(T+1,kneed)
        zlow = azlb
 347      pvmax = -0.1 d -9
*         JNEED = .18 d 0 / diffo
          JNEED = 10000
          xv(0) = 0.0 d 0
         DO 2010 jt = 1,JNEED
          xv(jt) = xv(jt-1) + diffo
         pval(jt) = pdfzp(xv(jt))
            tpv = pvmax - .1 d -4
            if (pval(jt).gt.tpv) then
               oldmax = pvmax
               pvmax = pval(jt)
               jmax = jt
               xmxa = pvmax * 0.0005 d 0
             else
               if (pval(jt).lt.xmxa) go to 1011
             end if
 2010  CONTINUE
 1011  trylas = xv(jt)
       table1(iptsig) = trylas
       table2(iptsig) = pval(jt)
         DO 102 iu = 1,jt
            if (pval(iu).gt.xmxa) go to 1021
  102    CONTINUE
 1021  trystr = xv(iu-1)
       table2(1) = pval(iu-1)
       iptsig = 50
       iepts = 78
       difg = (trylas - trystr) / dfloat(iptsig + iepts)
 8000  do 9907 itp = 1,iptsig
 9907 udifg(itp) = difg
       table1(1) = trystr
       table2(1) = 0.0 d 0
          udifg(2) = 20.0 d 0 * difg
          udifg(3) = 10.0 d 0 * difg
          udifg(4) = 6.0 d 0 * difg
          udifg(5) = 3.0 d 0 * difg
          udifg(49) = udifg(2)
          udifg(48) = udifg(3)
          udifg(47) = udifg(4)
          udifg(46) = udifg(5)
       do 276 ig = 2,iptsig
       table1(ig) = table1(ig-1) + udifg(ig)
 276   table2(ig) = pdfzp(table1(ig))
       prob = drnunf()

***** 
        do 7773, iv = 1, iptsig - 1
         if(table1(iv+1).le.table1(iv)) then
           write(50, 7990)
           write(50, 7991) kneed
           WRITE (50,1866) (theta(j),j=1,k),(sigf(jg), jg= 1,k),
     .  (kappa(jk), jk=1,K), (LAMBDA(JL),jl=1,K)
      WRITE (50,2066) (SIGMA(im,mj), mj = 1,M)
         do 2081 it = 0,T+1
 2081 write (50,2267) it,(z(it,nfa), nfa=1,k)
         go to 6646
           end if
 7773   continue
 6646   continue
*****

        z(T+1,kneed) = DGCIN(prob,2,iptsig,table1,table2)
 1      CONTINUE
        RETURN
        END
      DOUBLE PRECISION FUNCTION PDFZP(x)
*  subroutine to evaluate the kernel of z(T+1,kneed)
*    (For the 1 factor model, kneed is always 1.)
        IMPLICIT REAL * 8 (A-H,L,O-Z)
        DOUBLE PRECISION kappa,BSI(1000)
        INTEGER T
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /USE/ TWO,TOL
      COMMON /USFPAR/ A(500,5,3),B(500,5,3)
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /SOME/ SUM
         pi = 3.141592653589793238462643 d 0
* the only conditional: f(z(T+1) | z(T))
        vf = sigf(kneed) * sigf(kneed)
        dt = 1.0 d 0 - DEXP(-KAPPA(KNEED)*(1.0 d 0/52.14 d 0) )
        CIRC = (two * kappa(kneed)) / (vf * dt)
        CIRU = CIRC * z(T,kneed) * DEXP(-KAPPA(kneed) *
     e         (1.0 d 0/52.14 d 0) )
        CIRV = CIRC * x
        CIRQ = ( (two*kappa(kneed)*theta(kneed))/vf ) - 1.0 d 0
       ACIRQ = dabs(CIRQ)
      bbx = (CIRV / CIRU) ** (.5 d 0 * CIRQ)
      bx = two * ( (CIRU * CIRV)**.5 d0)
*     print 8650, bx
 8650 format (' bx: ',f12.5)
       if (bx.gt.1450.0 d 0) then
           zec = (8.0 d 0 * bx) * (8.0 d 0 * bx) * (8.0 d 0 * bx)
           BEST = 1.0 d 0 / dsqrt(two*pi*bx)
           bxmu = 4.0 d 0 * acirq*acirq
           best = best * (1.0 d 0 - (bxmu-1.0 d 0)/(8.0 d 0 * bx) +
     e ((bxmu-1.0 d 0)*(bxmu-9.0 d 0))/(two*(8.0 d 0*bx)*(8.0 d 0*bx))
     e - ((bxmu-1.0 d 0)*(bxmu-9.0 d 0)*(bxmu-25.0 d 0))/
     e  (6.0 d 0 *zec))
*         print 6785, acirq,best
 6785 format (' acirq: ',f12.3,' best: ',f12.4)
          best = dmax1(best,0.1 d -9)
      else
      if (ACIRQ.lt.1.0 d 0) then
         iu = 1
         go to 313
       end if
*     do 310 i = 2,200
*        if (dfloat(i).gt.ACIRQ.and.dfloat(i-1).lt.ACIRQ ) go to 312
*310   CONTINUE
*     Print 3166, ACIRQ
*3166 format (' problem - ACIRQ = ',f12.5)
* 312  iu = i-1
*      ACIRQ = ACIRQ - dfloat(iu)
       diu = dint(ACIRQ)
       iu = int(diu) + 1
         ACIRQ = DMOD(ACIRQ,diu)
 313   CALL DBSIES(ACIRQ,bx,iu,bsi)
* p. 10 in the special functions manual.
      BEST = BSI(iu)
      end if
 4565 bby = CIRC * DEXP(-CIRU-CIRV+bx)
      pdfzp = bby * bbx * BEST 
      if (pdfzp.lt.0.0 d 0) then
*        print 6579, bby,bbx,best,sum
 6579 format (' neggie pdfzP!',/' bby: ',f22.4,' bbx: ',f22.4,' best:
     f  ',f22.4,' sum: ',f22.4)
*        print 6678, acirq,bx,bxmu
 6678 format (' acirq: ',f22.6,' bx: ',f22.5,' bxmu: ',f22.4)
         STOP
       END IF
      RETURN
      END
      DOUBLE PRECISION FUNCTION PDFZ0(x)
*  subroutine to evaluate the kernel of z(0,kneed)
*   ineed is the point in time, kneed is the number of factor.
*    (For the 1 factor model, kneed is always 1.)
        IMPLICIT REAL * 8 (A-H,L,O-Z)
        DOUBLE PRECISION kappa,BSI(1000)
        INTEGER T
      COMMON /DATA/ R(500,5),Z(0:501,3),SMT(500,5),T,M,K
      COMMON /PARMS/ SIGMA(5,5),theta(3),kappa(3),LAMBDA(3),GAMMA(3),
     c sigf(3),SIGINV(5,5)
      COMMON /USE/ TWO,TOL
      COMMON /USFPAR/ A(500,5,3),B(500,5,3)
      COMMON /WHERE/ kneed,ineed,ngdr
      COMMON /SOME/ SUM
         pi = 3.141592653589793238462643 d 0
* The only conditional here: p(z(t+1) | z(t) )
        vf = sigf(kneed) * sigf(kneed)
        dt = 1.0 d 0 - DEXP(-KAPPA(KNEED)*(1.0 d 0/52.14 d 0) )
        CIRC = (two * kappa(kneed)) / (vf * dt)
        CIRQ = ( (two*kappa(kneed)*theta(kneed))/vf ) - 1.0 d 0
       ACIRQ = dabs(CIRQ)
        CIRU = CIRC * x * DEXP(-KAPPA(kneed) *(1.0 d 0/52.14 d 0) )
        CIRV = CIRC * z(1,kneed)
      bbx = (CIRV / CIRU) ** (.5 d 0 * CIRQ)
      bx = two * ( (CIRU * CIRV)**.5 d0)
       if (bx.gt.1450.0 d 0) then
           zec = (8.0 d 0 * bx) * (8.0 d 0 * bx) * (8.0 d 0 * bx)
           BEST = 1.0 d 0 / dsqrt(two*pi*bx)
           bxmu = 4.0 d 0 * acirq*acirq
           best = best * (1.0 d 0 - (bxmu-1.0 d 0)/(8.0 d 0 * bx) +
     e ((bxmu-1.0 d 0)*(bxmu-9.0 d 0))/(two*(8.0 d 0*bx)*(8.0 d 0*bx))
     e - ((bxmu-1.0 d 0)*(bxmu-9.0 d 0)*(bxmu-25.0 d 0))/
     e  (6.0 d 0 *zec))
           best = dmax1(best,.1 d -9)
      else
 313  if (ACIRQ.lt.1.0 d 0) then
         iu = 1
         go to 2313
       end if
*     do 2310 i = 2,200
*        if (dfloat(i).gt.ACIRQ.and.dfloat(i-1).lt.ACIRQ ) go to 2312
*2310   CONTINUE
*     Print 3166, ACIRQ
*3166 format (' problem - ACIRQ = ',f12.5)
*2312  iu = i-1
*      ACIRQ = ACIRQ - dfloat(iu)
       diu = dint(ACIRQ)
       iu = int(diu) + 1
       ACIRQ = DMOD(ACIRQ,diu)
 2313  CALL DBSIES(ACIRQ,bx,iu,bsi)
      BEST = BSI(iu)
* p. 10 in the special functions manual.
      END IF
*      print 1626, acirq,bx
 1626 format (' Bessel inputs: v: ',f10.7,' x: ',f22.8)
       cris = -ciru - cirv + bx
*      print *,' cris = ' , cris
 4565 bby1 = CIRC * DEXP(-CIRU-CIRV)
      bby = CIRC * DEXP(cris)
      pdfz0 = (bby * bbx * BEST)
*     print 35, x, pdfz0
*     print 36, bby,bbx,BEST
  35  format (' x: ',f12.7,' pdfz0: ',f30.14)
  36  format (' bby: ',f12.7,' bbx: ',f15.6,' best: ',f12.7)
      RETURN
      END